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Preface 


This  thesis  is  a  part  of  an  on-going  study  of  the  modelling  of  man-made  plasmas  under  Dr. 
William  F.  Bailey  at  the  Air  Force  Institute  of  Technology.  It  is  only  a  small  part  of  the  larger 
effort  to  which  Jones  [19],  Hilbun  [13],  Seger  [30],  Bell  [2],  Minor  [27],  and  many  others  have  made 
important  contributions. 

Because  this  work  is  part  of  a  larger  on-going  study,  I  have  included  in  it  two  items  that  are 
not  essential  to  this  project  itself,  but  should  be  of  help  as  a  reference  to  those  who  will  come  along 
later.  First,  the  detailed  derivation  of  the  moment  equations  in  Appendix  A  has  been  included 
to  provide  a  ready  reference  to  the  topic.  This  appendix  began  life  as  an  eight-page,  handwritten 
summary  in  Dr.  Bailey’s  introductory  plasma  physics  course  (which  in  turn  was  based  on  a  similar 
derivation  in  Holt  and  Haskell  [16:151-176]),  and  I  hope  that  it  will  make  the  full  circle  back  to 
being  a  class  handout;  it  answers  most  of  the  questions  that  I  had  as  I  was  trying  to  learn  the 
material.  The  second  example  of  extra  detail  is  the  taxonomy  of  reactions  in  a  hydrogen  ion  source 
in  Appendix  B.  Although  the  list  there  is  not  new  in  content,  I  hope  that  it  will  be  handier  in  form 
than  what  I  have  found  in  most  other  reports,  although  I  do  heartily  recommend  Janev  ei  al.  [17] 
as  a  complete  and  thorough  source  for  reaction  data.  Above  all,  I  hope  that  others  after  me  will  be 
able  to  find  some  foundation  here  upon  which  to  build,  as  I  have  built  upon  the  foundation  laid  by 
others  before  me.  To  this  end,  more  information  about  the  model  [35],  not  included  in  this  thesis, 
is  available  from  the  physics  department. 

A  special  thanks  goes  of  course  to  Dr.  William  F.  Bailey,  who  has  been  a  constant  source  of 
help  and  inspiration  throughout  this  project.  It  would  be  impractical  to  make  a  bibliographical 
reference: 

Bailey,  W.  F.  Personal  Interview. 
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because  it  would  appear  almost  everywhere  in  this  work,  so  I  gratefully  acknowledge  his  ubiquitous 
presence  here.  I  especially  appreciate  his  patience  with  a  student  who  probably  should  have  been 
a  music  major  during  a  term  that  seemed  much  more  like  a  race  against  the  clock  than  a  time 
for  thoughtful  pondering.  This  has  been  a  challenging  project  for  me,  which  has  included  learning 
from  scratch  FORTRAN,  UTeX,  emacs,  UNIX,  ftp,  etc. — not  to  mention  learning  a  lot  more  about 
plasma  physics  in  icn  sources.  Without  Dr.  Bailey’s  willingness  to  go  back  to  basic  principles  time 
and  time  again  when  I  didn’t  understand  what  was  going  on,  it  would  have  been  impossible. 

Finally,  I  owe  another  special  thanks  to  my  new  wife,  Ruth.  Since  our  marriage  she  has 
known  nothing  other  than  AFIT  widowhood,  yet  she  has  always  remained  patient  and  supportive, 
a  help  and  an  inspiration.  With  the  thanks  comes  a  promise  that  things  will  get  better  soon. 


This  thesis  was  typeset  in  DTeX  resident  on  AFIT’s  UNIX  VAX  computer  129.92.1.3,  nick¬ 
named  “Galaxy,”  The  model,  the  positive  ion  source  code  pos,  tan  on  Galaxy  in  UNIX  Fortran. 
The  MeiaLib  library  of  FORTRAN  routines  resident  on  Galaxy  provided  the  basis  of  the  graphics 
modules  that  graphed  the  results  of  the  model.  The  Macintosh  program  SuperPaint  drew  the  orig¬ 
inal  figures,  and  the  Macintosh  program  CrickeiGraph  redrew  the  experimental  data  plots  from 
Johnson  [18]. 


Todd  Roland  Vitk'^ 
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Abstract 

A  one-dimensional  fluid  model  of  pleisma  transport  in  tandem  volume  magnetic  multicusp  ion 
sources  is  explored.  The  model,  the  positive  ion  source  code  pos,  by  Glasser  and  Smith,  calculates 
plasma  density,  drift  velocity,  electron  temperature,  and  ion  temperature  in  an  ion  source.  The 
usefulness  of  the  model  is  limited:  (1)  The  plasma  density  trend  runs  opposite  to  experimental 
results,  and  electron  temperatures  are  an  order  of  magnitude  higher  than  experimentally  observed. 
(2)  Simplification  of  the  reaction  chemistry  leads  to  a  plasma  balance  between  ionization  and 
outflow  instead  of  the  correct  balance  between  ionization  and  recombination.  (3)  Wail  losses  are 
neglected.  (4)  There  are  inconsistencies  in  the  derivations  of  some  equations.  (5)  The  final  solution 
depends  on  the  choice  of  the  initial  estimated  solution.  (6)  Results  of  the  model  are  not  totally 
reproducible.  (7)  Numerical  instabilities  develop  upon  modification  of  terms  or  variation  of  initial 
conditions  outside  of  a  narrow  range.  Calculations  of  the  plasma  potential  from  the  results  of  the 
model  are  qualitatively  correct. 


NUMERICAL  ANALYSIS  OF  PLASMA 
TRANSPORT  IN  TANDEM  VOLUME 
MAGNETIC  MULTICUSP 
ION  SOURCES 


/.  Introduction 


The  purpose  of  a  neutral  particle  beam  is  to  transfer  energy  from  the  source  device  onto  a 
target.  One  of  the  primary  proposed  uses  of  neutral  particle  beams  is  use  as  a  weapon,  in  which 
case  the  transfer  of  energy  is  designed  to  destroy  the  target,  whether  it  be  an  incoming  nuclear 
warhead  or  a  hostile  space  platform.  The  other  principal  proposed  use  of  neutral  particle  beams 
is  for  the  injection  of  high-energy  neutral  particles  into  a  tokamakL  where  the  transfer  of  energy 
is  designed  to  heat  the  plasma  inside  the  tokamak  to  achieve  thermonuclear  fusion.  Since  interest 
in  the  use  of  particle  beams  as  weapons  has  waned  in  favor  of  other  options,  this  introduction  will 
motivate  the  present  study  in  terms  of  the  fusion  applications. 

For  a  hydrogen  plasma  to  achieve  fusion,  it  must  satisfy  the  Lawson  condition: 


nr  >  10‘® 

nr  >  10*“* 


d  —  d  reactions 
d  —  t  reactions 


}  T  >  IQkeV 


(1) 


with  sufficiently  high  density  n,  confinement  time  r,  and  temperature  T  [21:265].  One  way  to  heat 
the  plasma  to  achieve  the  Lawson  condition  is  through  radiofrequency  radiation;  another  way  is  to 
inject  particles  at  high  velocity.  The  high  directed  velocity  of  the  injected  particles  of  species  a 

'From  the  Russian  tokamak,  eui  acronym  for  toroidal’naya  kamrra  s  magnitnym  polem  “torroidal  cliamber  with 
a  magnetic  field.”  [28] 


1 


Q  is  transformed  (thermalized)  to  a  random  (thermal)  velocity  wf  through  collisions  with  particles 
of  the  plaisma.  To  make  injection  most  efficient,  a  beam  of  neutral  particles  is  injected. 

Injecting  neutral  particles  maintains  the  quasineutrality  of  the  plasma  and  thus  avoids  driving 
thermal  ions  to  the  wall.  However,  because  neutral  particles  are  uncharged,  they  can’t  be  acceler¬ 
ated  by  electric  or  magnetic  fields.  The  solution  to  this  impasse  is  to  accelerate  charged  particles 
and  then  convert  the  accelerated  particles  into  neutrals  without  significantly  reducing  their  average 
velocity. 

There  is  a  choice  between  two  types  of  charged  particles  to  accelerate:  positive  ions  or  negative 
ions.  Negative  ions  are  the  better  choice  for  three  reasons. 

One  advantage  that  negative  ion  beams  have  over  positive  ion  beams  is  that  the  negative  ion 
beams  are  more  easily  made  uniform.  For  fusion  applications,  the  injected  beam  should  ideally  be 
uniform.  Without  using  a  scheme  to  sort  by  mass,  any  beam  of  positive  ions  will  contain  the  species 
H2  and  H3  in  addition  to  //+.  Each  species  of  positive  ion  has  the  same  charge  and  therefore 
experiences  the  same  force  in  a  given  electromagnetic  field.  However,  because  each  species  has  a 
different  mass,  each  will  accelerate  at  a  different  rate  and  take  a  different  length  of  time  to  traverse 
a  given  distance.  So  any  unsorted  beam  of  positive  ions  will  have  three  components.  On  the  other 
hand,  ions  are  unstable  and  spontaneously  decay  into  their  constituent  components  [18:35],  so 
a  negative  ion  beam  contains  only  H~  ions  and  is  thus  uniform. 

Another  advantage  of  negative  ion  beams  over  positive  ion  beams  is  that  it  is  easier  to  turn  a 
negative  ion  beam  into  a  stream  of  neutral  particles  than  it  is  to  transform  a  positive  ion  beam  in 
the  same  way.  Since  in  the  negative  ion  beam  the  electrons  are  only  loosely  bound  to  the  hydrogen 
atoms,  they  are  relatively  easy  to  strip  off  of  the  atoms.  By  contrast,  it  is  relatively  difficult  to 
get  a  fast-moving  positive  ion  to  capture  an  electron.  Moreover,  the  cross  section  for  the  reaction 
//+  -\-e  ^  H  decreases  with  increasing  energy  [19:1-2],  so  the  more  energy  desired  to  be  input 
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into  the  tokamak,  the  more  difficult  it  will  be  to  neutralize  the  positive  ion  beam  that  carries  the 
energy. 

A  third  advantage  of  negative  ion  beams  over  positive  ion  beams  is  that  negative  ion  beams 
propagate  more  readily  through  the  drift  space  between  the  negative  ion  source  and  the  user  device. 
To  see  why  this  is  so,  consider  a  beam  of  negative  ions  propagating  through  a  low  density  gas.  Some 
of  the  ions  will  collide  with  the  atoms  of  the  gas,  producing  positive  ions  and  free  electrons.  Because 
the  electrons  have  a  much  lower  mass  than  the  positive  ions,  the  random  velocities  of  the  electrons 
will  be  higher  than  those  of  the  ions,  and  the  electrons  will  tend  to  leave  the  path  of  the  negative 
ion  beam  sooner.  The  departing  electrons  leave  behind  the  positive  ions,  which  form  a  sheath  of 
positive  charge  around  the  negative  ion  beam.  This  sheath  of  positive  charge  serves  as  a  conduit 
for  the  negative  ion  beam.  Thus  the  negative  ion  beam  is  said  to  be  “self-focusing”  [18:19]. 

Figure  1,  from  Johnson  [18:20],  shows  a  schematic  diagram  of  how  a  negative  ion  beam  could 
be  used  to  heat  a  tokamak.  Extraction  electrodes  extract  the  negative  ions  formed  in  the  source 
and  accelerate  them  towards  an  accelerator,  which  brings  them  to  their  maximum  velocity.  Before 
the  negative  ions  enter  the  tokamak,  they  peiss  through  a  plasma  neutralizer  which  strips  the  extra 
electrons,  transforming  the  beam  of  negative  ions  into  a  beam  of  neutral  particles  [H  atoms).  A 
magnetic  field,  located  just  past  the  neutralizer,  deflects  any  negative  ions  remaining  in  the  beam 
into  an  ion  dump. 

1.1  Negative  Ion  Sources 

The  majority  of  negative  ion  sources  may  be  classified  into  three  categories:  vapor  sources, 
surface  sources,  and  volume  sources. 

1.1.1  Vapor  Sources.  A  vapor  source  transforms  a  beam  of  positive  ions  into  a  beam  of 
negative  ions  by  sending  them  through  a  vapor  of  alkali  metal  atoms.  Because  the  alkali  metals 
have  an  electron  weakly  bound  in  the  outermost  shell,  this  electron  is  easily  stripped.  After  two 
captures,  a  positive  ion  in  the  incident  beam  is  transformed  into  a  negative  ion  in  the  exiting  beam. 
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Figure  1.  Neutral  Beam  Injection  (from  Johnson  [18:20]). 
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Figure  2.  Surface  Source,  after  Leung  and  Ehlers  [22:804]. 


However,  there  are  several  problems  with  this  kind  of  negative  ion  source.  First,  after  the 
incident  ion  captures  its  first  electron  and  before  it  captures  its  second,  it  is  a  neutral  atom.  There 
will  be  a  space  in  the  vapor  where  the  traversing  beam  will  be  made  up  mostly  of  neutral  particles. 
This  neutral  beam  will  diverge  before  the  second  capture  turns  the  neutral  atoms  into  negative  ions 
[18:21-22].  A  second  problem  with  the  vapor  source  is  high-voltage  breakdown  [18:22].  The  beam 
of  positive  ions  is  by  nature  an  electrical  current,  and  the  alkali  vapor  is  a  metal.  If  the  currents  are 
too  high,  the  metal  vapor  can  “short  out”  the  current  and  production  will  break  down.  The  third 
problem  is  that  collisions  with  alkali  atoms  will  impart  momentum  so  that  some  of  them  are  carried 
along  with  the  exiting  negative  particle  beam.  These  impurities  can  contaminate  the  plEisma  in  the 
tokamak  [18:22]. 


1.1.2  Surface  Sources.  A  second  source  of  negative  ions  is  the  cesiated  surface  source.  In 
this  source,  negative  ions  are  created  on  the  surface  of  a  cesiated  converter  electrode.  Figure  2 
shows  a  simplified  schematic  of  the  classical  surface  source  of  Leung  and  Ehlers  [22:804]. 
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In  the  surface  source,  electrons  emitted  from  hot  filaments  form  a  plaisma  in  the  vacuum 
vessel.  A  converter  electrode  in  the  center  of  the  vessel  maintains  a  large  negative  potential  (~  200 
V)  on  its  curved  front  face  and  a  small  positive  potential  (~  40  V)  on  its  sides  and  support  column. 
A  jet  injects  cesium  vapor  onto  the  surface  of  the  converter  electrode.  Attracted  by  the  negative 
converter  potential,  the  positive  ions  in  the  plasma  accelerate  toward  it  until  they  impact  the 
converter  surface.  Through  a  process  that  is  not  well  understood,  the  interaction  of  hydrogen  and 
cesium  ions  on  the  surface  produces  negative  hydrogen  ions.  These  newborn  negative  ions  are  now 
repelled  by  the  same  converter  potential  that  attracted  them  as  positive  ions.  Since  the  surface  of 
the  converter  electrode  is  curved,  the  ions  tend  to  accelerate  in  a  direction  perpendicular  to  the 
center  of  the  converter  electrode  and  are  thus  “focused”  geometrically  onto  the  extractor  electrode. 
Before  reaching  the  extractor  electrode,  the  ions  encounter  a  relatively  weak  (~  80  Gauss)  magnetic 
filter  field  which  returns  the  electrons  to  the  plasma  while  deflecting  the  negative  ions  only  slightly. 

There  are  also  several  problems  associated  with  the  surface  converter  source.  First,  because 
the  velocities  imparted  to  the  negative  ions  accelerating  away  from  the  converter  electrode  vary 
significantly,  the  beam  of  negative  ions  is  strongly  divergent  [18:22].  Second,  running  the  con¬ 
verter  electrode  potential  requires  an  additional  supply  of  power  [18:22].  Third,  the  relatively  large 
amounts  of  cesium  deposited  onto  the  converter  surface  and  onto  the  walls  of  the  vessel  create  a 
hazard  to  the  extraction  and  acceleration  system  [12:384]. 

The  primary  reaction  in  the  surface  source  depends  on  a  synergistic  interaction  between  the 
cesium  and  the  hydrogen.  The  H~  production  by  bombardment  with  both  cesium  and  hydrogen 
is  much  greater  than  the  production  due  to  either  process  alone  [31:432].  Yet  even  though  the 
primary  process  is  at  the  surface  of  the  converter  electrode,  no  less  than  25%  of  the  H~  ions  are 
produced,  not  at  the  surface  of  the  converter  electrode,  but  in  the  volume  of  the  device  [18:23]. 
With  such  a  high  production  in  the  volume  of  a  device  that  is  not  optimized  for  volume  production, 
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Figure  3.  Schematic  Diagram  of  a  Tandem  Volume  Source,  after  Glasser  and  Smith  [9:411]. 

the  natural  inclination  would  be  to  dispense  with  the  surface  converter  and  optimize  the  volume 
production  itself,  thus  creating  a  volume  source. 

1.1.3  Volume  Sources.  The  source  of  interest  to  this  study  is  the  tandem  volume  source,  or 
multipole  source,  depicted  in  Figure  3.  A  volume  source  generally  consists  of  a  cylindrical  vacuum 
vessel,  of  approximately  20  to  30  cm  in  length  and  20  to  30  cm  in  diameter,  divided  into  two  tandem 
chambers  by  a  magnetic  filter  field.  A  multipolar  array  of  permanent  magnets  surround  the  vessel 
to  conceal  the  walls  of  the  vessel  from  the  plasma,  limiting  losses  of  plasma  to  the  wall  to  the  cusps 
between  magnets  [20].  These  magnets  typically  have  strengths  of  about  3.6  kilogauss  (kG)  [23:56], 
The  inner  chamber,  called  the  source,  production,  or  driver  region,  contains  hot  filaments  that  emit 
fast  primary  electrons  (the  discharge).  These  filaments  form  the  cathode  of  the  device,  and  the  walls 
form  the  anode.  Typical  discharge  parameters  are  ~  80  volts  potential  and  ~10  amps  of  current 
[8:1425],  but  the  JET-type  positive  ion  source  at  the  Culham  laboratory  in  England  uses  discharge 
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currents  of  ~  1000  amps  to  optimize  H~  production  [18:24].  The  fast  primary  electrons  from  the 
discliarge  ionize  the  neutral  gas  in  the  source  chamber,  forming  a  plasma.  The  magnetic  filter  field 
confines  the  faist  electrons  to  the  source  region  but  allows  cold  secondary  electrons  to  diffuse  to  the 
second  chamber  [8:1423],  which  is  called  the  extraction  region  because  extraction  electrodes  at  the 
far  end  can  extract  ions  from  it.  The  volume  source  can  produce  either  negative  ions  or  positive 
ions,  depending  on  the  gas  mixture  used  and  whether  the  extraction  electrode  supplies  a  positive 
or  negative  potential  relative  to  the  plasma  potential. 

There  are  two  means  of  generating  the  magnetic  filter  field  [18:61].  Figure  4  shows  a  photo¬ 
graph  of  the  tandem  volume  source  at  Lawrence  Berkeley  Laboratory,  taken  from  Stevens  [33:272]. 
This  source  is  an  example  of  the  “rod-filter”  configuration,  where  magnetic  rods  physically  intrude 
across  the  interior  of  the  source  to  generate  the  magnetic  filter  field.  The  other  means  of  generating 
the  filter  field  is  by  rearranging  the  multipolar  array  of  magnets  on  the  outside  of  the  vessel  to 
extend  a  pair  of  cusps  so  that  they  cross  the  volume,  creating  a  “virtual”  magnetic  filter  field. 
Figure  5,  from  McAdams  [26:303],  shows  a  line  drawing  of  a  tandem  volume  source  with  a  virtual 
magnetic  filter. 

Chapter  II  discusses  the  kinematics  of  tandem  volume  ion  sources  in  greater  detail. 

J.2  Problem  Statement  and  Objectives 

The  purpose  of  this  thesis  is  to  investigate  the  utility  of  a  one-dimensional  plasma  fluid  model 
in  the  study  of  tandem  volume  magnetic  multicusp  ion  sources.  The  objectives  were:  (1)  Investigate 
the  kinematics  of  hydrogen  ion  sources  and  apply  what  was  learned  to  understand  the  model  under 
study.  (2)  Determine  the  origin  and  the  validity  of  the  equations  on  which  the  model  is  based. 
(3)  Examine  the  results  of  the  model,  especially  the  effect  of  the  magnetic  filter  field.  (4)  Compare 
the  results  from  the  model  to  what  would  be  expected  theoretically  and  experimentally  from  a 
positive  ion  source.  (5)  Where  the  results  differ  from  what  would  be  expected,  propose  a  reason 


8 


each  topic  (some  topics  are  more  heavily  weighted  to  one  category  than  another).  This  intro¬ 
ductory  chapter  concludes  with  a  short  presentation  of  the  salient  features  of  the  model  and  an 
explanation  of  the  notation  used  in  this  thesis.  Chapter  II  contains  a  brief  general  and  theoretical 
exposition  of  the  kinematics  of  hydrogen  ion  sources.  Chapter  III  applies  the  kinematics  discussed 
in  Chapter  II  to  evaluate  the  plasma  balance  in  the  model.  Chapter  IV  examines  the  results  of 
the  model’s  published  '‘test  case”  [9:422-423]  and  identifies  behavior  not  mentioned  in  the  original 
article.  Chapter  V  attempts  to  explain  some  of  the  behavior  discussed  in  Chapter  IV'  in  terms  of 
the  numerical  processes  occuring  in  the  model.  Chapter  VI  summarizes  the  development  of  the 
equations  used  in  the  model.  Chapter  VII  identifies  some  apparent  inconsistencies  in  the  derivation 
outlined  in  Chapter  VI  and  presents  the  results  of  modifications  to  the  model  designed  to  remove 
the  inconsistencies.  Chapter  VIII  describes  four  enhancements  to  the  model.  The  chapter  for 
overall  conclusions  and  recommendations  for  further  study  comes  in  the  usual  position.  The  math¬ 
ematical  details  of  the  derivation  of  the  moment  equations  discussed  in  Chapter  VI  are  relegated 
to  Appendix  A.  Appendix  B  gives  one  possible  taxonomy  of  reactions  in  a  hydrogen  plasma,  based 
on  the  presentation  in  Smith  and  Glasser  [32:393-394],  Appendix  C  is  a  glossary  of  symbols. 

1.4  The  Model 

The  model  used  in  this  investigation  is  the  positive  ion  source  code  pos,  developed  in  1987  by 
Alan  H.  Glcisser  and  Kenneth  Smith  at  the  Los  Alamos  National  Laboratory,  which  applies  fluid 
equations  to  solve  for  the  time  evolution  of  plasma  density,  drift  velocity,  electron  temperature, 
and  ion  temperature  in  a  positive  ion  source.  It  is  one  step  in  a  program  to  develop  a  “realistic 
computer  model  of  H~  ion  sources”  [9:410].  Although  the  physical  setup  of  the  negative  ion  source 
is  identical  to  that  of  the  positive  ion  source,  the  number  of  reacting  species  and  the  equation  set 
that  the  model  uses  to  describe  the  positive  ion  source  is  smaller  than  that  required  to  model  a 
negative  ion  source.  Specifically,  the  model  considers  only  four  interacting  species:  e,  ,  H .  and 
//o.  Thus  it  does  not  model  a  negative  ion  source  because  it  excludes  H~  ions  altogether. 
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The  pos  code  can  be  described  £is  a  physics  “chassis”  built  around  the  numerical  engine  dirk2. 
The  dirk2  code,  written  by  Neil  Carlson  and  Keith  Miller  of  the  University  of  California  at  Berkeley 
in  1984,  solves  a  system  of  first-order,  linearly-implicit,  ordinary  differential  equations  of  the  form 

A{y)  ^  =  9{y) 

using  a  diagonally  implicit  ^unge-J^utta  method  of  _^d  order  and  Miller’s  method  of  moving  finite 
elements  [5].  The  pos  code  establishes  the  physics  and  casts  the  differential  equations  in  a  form 
that  dirk2  can  solve. 

As  FORTRAN  code,  the  pos  code  is  well-written,  well-documented,  and  portable.  Although 
the  code  was  written  for  a  Cray  supercomputer,  it  compiled  without  error  in  UNIX  FORTRAN  on 
the  UNIX  VAX  “Gala.xy.”  Even  the  extremely  sensitive  VMS  FORTRAN  compiler  on  the  VMS 
VAX  “Hercules”  detected  only  one  small  error  during  compilation:  the  local  counter  variable  treac 
in  the  module  c/irate  has  the  same  name  as  the  array  ireocfnreocj  in  the  common  block  fesrofe.  This 
common  block  is  available  to  the  chraie  module,  but  the  array  ireac(nreac)  is  not  used  in  chrate. 

Figure  6  shows  a  block  diagram  of  the  subroutines  of  the  pos  code.  These  subroutines  may 
be  divided  into  two  groups.  The  first  group  consists  of  those  subroutines  which  are  either  provided 
with  the  dirk2  code  or  required  by  it.  The  dirk2  code  provides  the  routines  dirk2,  drkpfl,  bee,  and 
tnirp,  and  it  requires  the  user  to  supply  the  routines  resdl,  jac,  yprime,  linsoi,  norm,  setbvs,  solchk, 
and  step.  Also  numbered  among  this  first  group  are  the  subroutines  vdec,  vsol,  scopy,  solbt,  deebi, 
amaxaf,  display,  timer,  and  exit,  which  are  required  by  norm,  linsol,  and  jac.  The  second  group  of 
subroutines  are  those  which  establish  the  physics  of  the  ion  source.  This  second  group  is  comprised 
of  those  subroutines  that  either  input  and  initialize  the  physical  parameters  {rates,  xsecn,  and  inif) 
or  actually  compute  the  physics  (arrays,  fluxes,  brag  chem,  chraie,  reloss,  and  crosec).  This  thesis 
deals  almost  exclusively  with  the  second  group  of  subroutines,  which  are  diagrammed  separately 
in  Figure  7. 
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Figure  6.  Block  Diagram  of  the  pos  Code,  Showing  Subroutine  Calls. 
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main 


Figure  7.  The  Heart  of  the  Physics  of  the  pos  Code. 

The  main  module  of  pos  reads  the  initial  conditions  stored  in  the  file  pos.in,  defines  the  initial 
parameters  required  by  dirkS,  initializes  and  scales  the  variables,  and  writes  the  output  to  the  file 
pos. out.  It  also  contains  a  loop  of  subroutine  calls  that  accomplish  the  actual  calculation. 

As  part  of  the  process  of  initialization,  pos  calls  the  module  rates,  a  subroutine  that  reads  the 
chemical  reaction  rate  data  stored  in  the  file  chrate.i  and  stores  it  into  the  common  block  ksrate,  and 
the  module  rsecn,  a  subroutine  that  reads  the  elastic  cross  section  data  stored  in  the  file  crosec.i 
and  stores  it  into  the  common  block  coefts.  These  two  input  files  begin  two  parallel  calculations, 
one  for  the  chemical  reaction  rates,  the  other  for  the  elastic  cross  sections,  which  converge  in  the 
module  fluxes. 

The  module  chrate  is  a  subroutine  that  uses  the  variables  stored  in  the  common  block  ksrate 
to  calculate  the  reaction  rates  for  the  chemical  reactions  modelled  by  pos.  It  stores  these  reaction 
rates  in  the  variable  array  rates  to  pass  to  the  module  chem. 
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The  parallel  module  crosec  is  a  subroutine  that  uses  the  variables  stored  in  the  common  block 
coefts  to  calculate  the  elastic  scattering  cross  sections  with  neutral  particles.  It  passes  these  cross 
sections  to  the  module  brag  in  the  variable  array  sigma. 

The  chem  module  uses  the  reaction  rates  it  receives  from  the  chrate  to  calculate  the  collisional 
source  terms  5,  ii*",  Q*",  and  Q*",  which  it  stores  in  positions  1  through  5  of  the  variable 

array  sc,  which  it  passes  to  fluxes. 

The  parallel  calculation  in  the  brag  module  uses  the  el^lstic  scattering  cross  sections  it  receives 
from  the  crosec  a  number  of  plasma  parameters  and  coefficients,  which  are  stored  in  arrays  in  the 
common  block  coef,  which  is  available  to  the  module  fluxes. 

The  two  parallel  streams  of  calculations  now  converge  in  the  module  fluxes.  The  subroutine 
fluxes  uses  the  collisional  source  terms  from  ckem  and  the  plasma  parameters  and  transport  coeffi¬ 
cients  from  brag  to  calculate  convective  and  dissipative  flux  and  source  terms  that  are  components 
of  the  fluid  equations. 

The  fluxes  module  feeds  these  terms  to  the  module  arrays,  which  is  a  subroutine  that  builds 
the  arrays  that  will  be  manipulated  by  the  matrix  arithmetic  modules  and  eventually  solved  by 
dirk2. 

All  of  the  subroutines  are  further  documented  in  the  source  code  [5],  in  the  original  article  [9], 
or  in  an  additional  report  [35]. 

1.5  Notation 

Notation  is  a  difficult  subject  in  plasma  physics.  There  are  too  few  letters  in  the  Latin  and 
Greek  alphabets  to  denote  the  many  different  quantities  that  can  be  defined.  Moreover,  notation 
is  not  quite  standard,  and  sources  differ  considerably.  This  thesis  attempts  to  fulfill  the  (not 
always  consistent)  twin  goals  of  (1)  using  notation  that  is  widely  recognized  and  (2)  using  the  same 
representation  throughout  for  any  given  quantity. 
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This  thesis  uses  indicial  notation  throughout.  On  the  few  occasions  where  a  vector  is  not 
written  with  indices,  it  is  put  into  bold  type:  H,  V“.  Instead  of  the  Greek  letters  a  and  /?  that 
Braginskii  uses  for  indices  [3:205],  this  thesis  uses  the  Latin  letters  i,  j,  and  ^  as  in  Holt  and  Haskell 
[16]  and  in  Golant  [10].  Each  index  represents  the  principal  coordinate  directions  x,  y,  and  z.  The 
ion  source  geometry,  depicted  in  Figure  3,  defines  these  principal  coordinate  directions.  Moving 
from  the  production  region  to  the  extraction  region  along  the  axis  of  the  ion  source  defines  the 
positive  z  direction,  with  zero  at  the  filament  end.  The  positive  direction  of  the  magnetic  filter 
field  defines  the  positive  direction  of  the  x  axis.  The  right-hand  rule  determines  the  positive  sense 
of  the  y  axis. 

Subscripts  are  reserved  for  (1)  vector  indices  and  (2)  integration  bin  numbers.  Thus  t;,-  denotes 
a  component  of  velocity  in  the  direction,  and  nj  the  plasma  density  in  the  integration  bin. 
The  context  makes  clear  which  convention  is  which. 

To  avoid  confusion  with  the  subscripts  for  indices  and  integration  bin  numbers,  species  desig¬ 
nations  are  superscripts.  In  the  two-fluid  derivations,  i  designates  ions  and  e  electrons.  Otherwise, 
species  are  designated  with  their  chemical  identifications,  c,  H~ ,  H,  H2,  ,  and  .  Thus 

is  the  density  of  hydrogen  atoms,  the  component  of  average  ion  velocity  in  the  z  direction, 
and  T'  the  electron  temperature.  Instead  of  the  Latin  letter  a  that  Braginskii  uses  [3:205]  or  the 
Latin  letter  j  that  Glasser  and  Smith  use  [9:412]  to  designate  a  general  (charged)  species,  this  thesis 
uses  the  Greek  letter  o  in  accordance  with  Krall  and  Trivelpiece  [21:79]  and  Golant  [10],  Thus 
is  the  density  of  particles  of  species  o.  As  in  Golant  [10:315],  the  Greek  letter  0  designates  a  second 
general  (charged)  species:  for  example,  the  general  momentum  exchange  term  ±.R°^  represents 
the  momentum  transferred  from  ions  to  electrons,  -t-i?**,  in  the  ion  equation  (where  q  represents 
species  i,  so  0  represents  the  overcharged  species  c),  but  it  represents  the  momentum  transferred 
from  electrons  to  ions  —R‘'  in  the  electron  equation  (where  the  general  charged  species  q  is  e,  and 
the  other  charged  species  0  becomes  i). 
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Unfortunately,  superscript  species  designations  look  confusingly  like  standard  notation  for 
exponents  (which  are  seldom  used  in  this  thesis).  To  mitigate  confusion  where  it  could  arise, 
exponents  are  written  outside  of  parentheses.  Thus  the  square  of  the  velocity  vector  would  be  (v)^. 
However,  in  keeping  strictly  with  indicial  notation,  this  thesis  uses  the  notation  r,u,-  instead  of  (v)* 
whenever  possible. 

Since  reaction  rates  and  cross  sections  are  never  subscripted  with  integration  bin  numbers 
or  vector  directions,  the  chemical  reaction  numbers  in  Appendix  B  are  written  as  subscripts.  The 
Hebrew  letter  N  (aleph)  denotes  the  general  reaction.  Thus  kis  designates  the  reaction  rate  for 
Reaction  19  and  (Th  the  cross  section  for  Reaction  K. 

For  clarity,  partial  derivatives  are  written  out  in  full,  not  abbreviated  with  commas.  Hence, 

dV°' 

the  time  derivative  of  the  mean  velocity  in  the  i‘''  direction  is  ,  not  . 

at 

Lastly,  Braginskii  denotes  a  cross  product  with  square  brackets,  [vB]o  [3:205].  This  thesis 
uses  the  indicial  convention  for  the  alternating  unit  tensor  defined  in  Appendix  A,  djkVjBk,  but  as 
a  concession  to  Braginskii,  all  cross  products  will  also  be  enclosed  in  brackets  for  clarity:  [tijkVjBk]- 

Appendix  C  summarizes  the  nomenclature  used  in  this  thesis. 
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II.  The  Kinematics  of  Negative  Ion  Sources 


In  the  hydrogen  plasma  in  the  ion  source,  there  are  seven  interacting  species,  which  are  listed 
in  Table  1. 


mass  — » 

charge 

1 

c  H- 

H  Hi 

H+  Hf  Hf 

Table  1.  The  Seven  Species  in  a  Hydrogen  Ion  Source. 

There  are  28  possible  pairwise  interactions  among  the  seven  species.  In  addition,  each  pairwise 
interaction  can  in  general  generate  several  different  possible  combinations  of  reaction  products — 
the  more  subatomic  components  in  the  reactants,  the  greater  the  possible  variety  of  products  of 
a  given  interaction.  Moreover,  all  of  the  reactant  species  except  H'^  and  e  have  internal  energy 
states  [34:15],  such  as  the  vibrationally  excited  states  (4  <  v"  <  11)  of  H2  or  the  electronically 
excited  state  H‘{2p)  of  H.  The  symbol  v"  denotes  a  vibrationally  excited  state,  and  the  symbol 
a*  denotes  an  electronically  excited  state  of  species  o.  In  calculating  the  reaction  chemistry,  most 
of  the  excited  states  must  be  treated  as  separate  particles  because  the  cross  sections,  and  thus  the 
chemical  reaction  rates,  differ  from  those  of  the  unexcited  particles  [34:17].  All  things  considered, 
the  taxonomy  of  reactions  in  even  the  simplest  hydrogen  plasma  is  quite  complex.  Appendix  B 
contains  basic  lists,  after  Smith  and  Glasser  [32:393-394],  of  the  principal  reactions  that  form  and 
destroy  the  seven  species  under  consideration. 

Although  tandem  magnetic  multicusp  ion  sources  were  originally  developed  as  positive  ion 
sources  [18:23],  with  the  magnetic  filter  field  serving  to  alter  the  species  ratio  ’Hf)  of  the 

positive  ions  extracted  [25:1],  currently  most  are  constructed  to  optimize  the  production  of  H~ 
ions. 

The  kinematics  of  the  tandem  ion  source  (and  hence  the  production  of  H~  ions)  are  deter¬ 
mined  to  a  large  extent  by  (1)  the  production  of  plasma  through  ionization  at  the  hot  filament 
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cathode,  (2)  the  losses  of  plasma  (a)  to  the  walls  of  the  chamber,  (b)  through  recombination  in  the 
volume,  and  (c)  out  the  extraction  end  of  the  device,  and  (3)  the  properties  of  diffusion  through 
the  magnetic  filter  field. 

2.J  Plasma  Production  at  the  Cathode 

To  provide  the  material  with  which  to  form  a  plasma,  a  gas  jet  admits  neutral  hydrogen  gas 
into  the  chamber.  A  cathode  consisting  of  hot  filaments  (~  100  eK)  supplies  energy  to  the  gas  in 
the  source  chamber  in  the  form  of  fast  electrons  emitted  from  the  filaments.  These  fast  primary 
electrons,  e^"^  ,  form  what  is  probably  an  offset  (“drifting”)  Maxwellian  distribution  [18:28]  of 
energies  and  velocities  around  the  injection  potential  of  the  filaments.  The  primary  electrons 
ionize  molecules  of  neutral  hydrogen  gas,  losing  energy  in  each  collision  until  they  cross  eth ,  the 
threshold  energy  for  the  cross  section  for  ionization,  where  they  join  a  Maxwellian  population  of 
thermal  secondary  electrons,  Because  the  fast  primary  electrons  give  up  their  energies  to  the 

discharge  and  join  the  population  of  secondary  electrons  much  more  rapidly  than  the  secondary 
electrons  are  lost  to  the  walls  or  are  extracted  out  of  the  device,  the  population  of  thermal  secondary 
electrons  grows  to  be  several  orders  of  magnitude  larger  than  the  population  of  primary  electrons. 

Ionization  is  not  the  only  process  that  redistributes  the  energy  of  the  primary  electrons.  Elec¬ 
tron  excitation  and  vibrational  excitation  also  redistribute  energy,  and  each  has  its  eissociated  cross 
section.  Figure  8  shows  a  schematic  diagram  of  the  interaction  of  the  electron  energy  distribu¬ 
tion  function  (EEDF)  with  the  various  cross  sections.  This  digram  shows  that  primary  electrons 
accomplish  most  of  the  ionization  of  the  neutral  gas  to  form  a  plasma. 

Primary  electrons  may  have  energies  up  to  the  discharge  potential  [18:28],  which  is  commonly 
of  the  order  of  lOOV  [4:379].  The  higher  the  energy  of  the  primary  electron,  the  more  electron-ion 
pairs  it  can  produce  before  it  joins  the  population  of  thermal  electrons  [24:365].  Suppose  that  a 
narrow  beam  of  electrons  impinges  on  a  plasma  so  that  the  incident  electrons  have  a  narrow  range 
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Figure  8.  Schematic  Diagram  of  the  Relative  Range  of  Overlap  of  Typical  Cross  Sections  and  the 
Electron  Energy  Distribution  Function. 

of  energies  centered  around  Cb«am.  as  depicted  in  Figure  9.  Then  an  electron  in  the  beam  will 
lose  energy  through  collisions,  mainly  with  neutral  particles.  Often  it  will  ionize  these  particles 
in  the  course  of  the  interaction.  The  greatest  possible  number  of  secondary  electrons  would  be 
produced  if  (1)  every  interaction  caused  an  ionization  and  (2)  each  of  the  secondary  electrons  thus 
produced  were  born  at  rest  (a,  6,  c,  d,  e,  f,  and  g  in  the  diagram).  In  this  case,  all  of  the  energy 
lost  by  the  primary  electron  goes  into  breaking  the  electron  bond  and  none  is  carried  away  by 
the  secondary  electron  as  kinetic  energy.  Thus  the  primary  electron  would  produce  the  greatest 
number  of  electron-ion  pairs  before  its  energy  falls  below  the  threshold  energy  Cth.  at  which  the 
collision  cross  section  effectively  goes  to  zero.  A  100  eV  primary  electron  could  accomplish  seven 
“generations”  of  13.58  eV  ionizations  before  falling  below  the  13.58  eV  threshold  energy.  On  the 
other  hand,  if  the  ionized  secondary  electrons  were  born  with  a  substantial  kinetic  energy  Ckinet 
{a,  0,  7,  and  S  in  the  diagram),  then  the  number  of  ionizations  that  the  primary  electron  could 
accomplish  before  joining  the  population  of  thermal  secondary  electrons  would  be  much  smaller  («  A 
for  Ckinet  »  10  eV).  The  truth  apparently  lies  somewhere  between  the  two  extremes:  modeling  by 
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Chan,  Burrell,  and  Cooper  shows  that  the  first  generation  of  primary  electrons  accomplish  roughly 
half  of  the  ionizations  [7:6124],  yet  it  takes  a  full  five  generations  to  convert  a  population  of  primary 
electrons  born  at  $  as  60  V  into  secondary  electrons  [7:6123], 


2.2  Wall  Losses 

Wall  losses  and  recombination  are  the  two  principal  losses  of  charged  particles  in  a  tandem 
volume  multicusp  ion  source  [7:6124],  As  ions  reach  the  wall,  which  also  serves  as  the  anode  of 
the  discharge,  they  are  generally  neutralized.  Although  the  wall  is  protected  with  a  short-range 
multipole  magnetic  field,  ions  are  lost  to  the  cusps  between  the  magnets.  The  apparent  width  of 
the  cusps  varies  according  to  species,  appearing  the  smallest  to  the  fast-moving  primary  electrons, 
rather  larger  to  the  slow  secondary  electrons,  and  larger  still  to  the  ions  [18:31],  Notwithstanding 
the  apparent  size  of  the  cusps,  the  higher  the  energy  of  the  primary  electrons,  the  higher  their  loss 
rate  to  the  walls  [24:365], 
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In  addition  to  neutralizing  electrons  and  H  ions,  the  walls  will  also  neutralize  H  ions. 
Johnson  [18:42]  estimates  the  rate  constant  ku,aU  for  neutralization  of  H~  in  the  reaction 

H~  -f  (wall)  — ►  /f  +  (wall) 

in  terms  of  the  plasma  potential  ^pia.ma  and  the  negative  ion  temperature  s«  0.4  eV  to  be 

Jb^aH  =  2.4x  (2) 

One  could  argue  that  a  one-dimensional  model  would  lose  plasma  to  the  wall  only  at  the  inner 
boundary  «  =  0.  However,  losses  at  the  wall  along  the  entire  length  of  a  real  ion  source  lower  the 
plasma  density  near  the  wall,  and  the  high  random  thermal  velocity  of  the  plasma  quickly  transmits 
the  effect  of  this  lowered  density  near  the  wall  throughout  the  volume.  Therefore,  a  rigorous  model 
must  account  for  wall  losses  in  some  fashion. 

2.3  The  Magnetic  Filter  Field 

The  magnetic  filter  field  divides  the  tandem  ion  source  into  driver  and  extraction  chambers. 
In  the  axial  direction,  its  profile  is  approximately  Gaussian,  with  an  intensity  almost  two  orders 
of  magnitude  lower  than  the  multipolar  field  protecting  the  walls  (~  80  Gauss  vs.  ~  3600  Gauss). 
The  V  X  B  force  that  the  charged  particles  experience  in  the  field  causes  them  to  orbit  around  a 
“guiding  center”  at  a  distance  equal  to  their  Larmor  radius,  ,  given  by 


“  I  I  c5 

Because  the  ions  have  much  greater  mass,  their  gyroradius  in  the  weak  magnetic  filter  field  is  so 
large  that  they  are  deflected  only  slightly  as  they  pass  through  it  [23:56].  Electrons,  on  the  other 
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hand,  are  either  deflected  by  the  filter  field  and  returned  to  the  extraction  region  discharge  (if  their 
energies  are  sufficiently  high)  or  are  trapped  by  the  magnetic  field  lines.  As  these  trapped  electrons 
orbit,  they  collide  with  other  particles  in  the  volume.  Each  collision  causes  the  guiding  center  of 
the  electron’s  orbit  to  be  displaced  randomly,  so  the  electron  begins  a  “random  walk”  diffusion 
in  the  magnetic  filter.  The  electrons  that  reach  the  extraction  chamber  through  diffusion  are  all 
thermalized  (cold)  from  collisions.  Thus  the  magnetic  filter  effectively  creates  a  barrier  that  stops 
hot  primary  electrons  but  allows  cold  secondary  electrons  to  pass  through  readily  [8:1423],  limiting 
the  mean  electron  temperature  in  the  extraction  chamber  to  about  0.4 eV^  [23:56].  A  magnetic 
filter  field  of  strength  B^30  Gauss  is  sufficient  to  prevent  the  primary  electrons  from  reaching  the 
extraction  region  [18:295-296]. 

2.^  The  Formation  of  Negative  Ions 

Negative  ions  are  formed  in  a  two-step  process.  In  the  first  step,  fast  primary  electrons 
excite  electronic  states  of  precursor  molecules,  which  spontaneously  relax  to  vibrationally  excited 
states.  In  the  second  step,  slow  electrons  produce  negative  ions  from  the  vibrationally  excited 
molecules  through  dissociative  attachment  reactions.  The  tandem  multicusp  volume  source,  wnth 
two  chambers  separated  by  a  magnetic  field,  is  ideally  suited  to  meet  these  inherently  incompatible 
production  requirements. 

The  primary  route  of  formation  of  negative  ions  in  a  volume  of  plasma  is  through  dissociative 
attachment  of  an  electron  to  a  vibrationally  excited  gas  molecule  [18:38]: 

e  + //2  (v") -*//"  + //  (3) 

where  v"  signifies  a  vibrationally  excited  state  of  H2-  This  reaction  is  labelled  Reaction  56  in 
Appendix  B.  Reaction  56  obviously  requires  a  large  population  of  H^iv")  to  create  a  significant 
amount  of  H~ .  This  is  where  the  tandem  volume  source  has  its  great  advantage.  Creating  sufficient 


amounts  of  H^W)  requires  a  population  of  fast  electrons;  however,  these  same  fast  electrons  that 
form  H2W)  neutralize  H~  ions  at  a  high  rate  by  way  of  the  collisional  detachment  reaction: 

e  -p  //  — *-  e  +  //  -f  e 

which  is  Reaction  7  in  Appendix  B.  In  the  tandem  volume  source,  the  driver  region  is  populated 
with  fast  primary  electrons  that  can  excite  the  vibrational  modes  of  H2,  and  the  extraction  region 
is  free  from  these  same  fast  electrons  that  would  dissociate  H~ .  According  to  Hiskes  [14:3],  there 
are  three  sources  of  (1)  relaxation  of  electronically  excited  by  collisions  with  primary 

electrons,  (2)  neutralization  of  /f^at  the  wall,  and  (3)  relaxation  of  electronically  excited  by 
collisions  with  the  wall. 

The  sequence  of  events  runs  roughly  as  follows:  In  the  driver  region,  a  primary  electron, 
gives  up  kinetic  energy  to  excite  an  H2  molecule  electronically  [15:34]; 

€^"^fast)  +  Hi  ^  +  H; 

The  electronically  excited  molecule  next  relcixes  to  a  vibrationally  excited  state  Hi{v")  [18:35] 

Approximately  90%  of  the  H~  ions  extracted  from  the  source  are  derived  from  Hi  {5  <  v"  <  11) 
[14:6],  Because  the  cross  sections  are  smaller  for  the  lower  vibrational  states,  and  because  the 
probability  of  an  electron  spontaneously  detaching  from  the  ion  increases  significantly  at  vibrational 
excitation  levels  above  v"  =  9  [18:36],  the  most  important  vibrational  states  are  6  <  r"  <  9  [18:36]. 
The  increase  in  probability  of  reaction  with  vibrational  states  is  significant:  the  cross  section  for 
Reaction  56  in  Equation  3  increases  about  four  orders  of  magnitude  from  v"  =  0  to  r"  =  4  [18:36]. 
The  excitation  may  actually  proceed  in  a  two-  or  three-step  process,  since  the  probability  for  the 
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H-i  molecule  to  reach  5  <  v"  <  11  through  collision  with  an  electron  is  much  greater  when  it  begins 
from  the  vibrational  population  at  v"  =  1  or  v"  =  2  instead  of  the  one  at  v’'  =  0  [15:34]. 

The  magnetic  filter  field  is  transparent  to  neutral  H2  molecules,  so  they  readily  reach  the 
extraction  region,  where  they  collide  with  slow  secondary  electrons  and  can  form  H~  by  way 
of  the  intermediate  unstable  molecular  ion  [18:35]; 

H  -f-  II~ 

(slow)  +  //a  (5  <  v"  <  11)  /f-  J 

^2  +  e 

In  the  extraction  region,  the  principal  loss  processes  for  H~  are  the  mutual  neutralization  Reaction 
11  with  ions  and  the  associative  detachment  Reaction  58  with  H  atoms  [14:6]: 

—  H  +  H  and  ^  + /f"  ^  f/a  (v")  +  e 

The  chain  of  reactions  that  produce  H~  appears  to  be  quite  successful  because  in  the  extraction 
region  the  density  of  H~  can  approximate  the  density  of  electrons  [18:35-36]. 

2.5  Summary  of  Kinemattcs 

Figure  10  summarizes  the  relevant  kinematic  processes  in  a  tandem  multicusp  volume  hydro¬ 
gen  ion  source.  A  jet  admits  hydrogen  gas  into  the  chamber,  and  a  cathode  inputs  fast  primary 
electrons  into  the  driver  region.  These  primary  electrons  ionize  the  hydrogen  gas,  forming  a  plasma. 
Recombination  reactions  among  the  reactants  and  with  the  wall  limit  the  plasma  density  and  pro¬ 
duce  vibrationally  excited  gas  molecules.  The  magnetic  filter  field  transmits  neutral  atoms  and 
molecules,  positive  ions,  and  cold  secondary  electrons,  while  acting  as  a  barrier  to  fast  primary 
electrons.  In  the  extraction  region,  interactions  between  cold  electrons  and  vibrationally  excited 
gas  molecules  create  negative  ions  while  other  interactions  with  the  reactants  tend  to  destroy  them. 


25 


Figure  10.  Schematic  of  the  Kinematics  of  a  Tandem  Volume  Hydrogen  Ion  Source,  from  Johnson 
[18:29], 
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III.  The  Plasma  Balance 


Glctsser  and  Smith  based  the  model  on  a  two-fluid  and  e)  derivation  of  the  plasma 

transport  equations,  enhanced  by  the  effect  of  collisions  with  neutral  species  {H  and  H2)  [9:409]. 
Thus  of  the  seven  species  in  Table  1  in  Chapter  II,  the  model  considers  only  four  (e,  H'^ ,  H ,  and 
H2)  and  reduces  the  number  of  reactions  to  seven,  as  given  in  Table  2. 


Inelastic  Chemical  Reactions 

ku  = 

N 

n 

reaction 

name 

£th 

TJ  =  10  eU 

19 

e 

e  +  H 

— ► 

H+  +  2e 

ionization 

5.7  X  lO-"' 

23 

i 

H  +  H+ 

— 

H+  +  H 

charge  exchange 

iB 

QtR  gn 

3.2  X  10'® 

26 

e 

e  +  H+ 

— 

H  +  hu 

radiative  capture 

0.0 

Ski 

3.2  X  10"^^ 

Elastic  Scattering  Reactions 

K 

j 

reaction 

name 

^th 

-leV 

34 

e 

e  -f-  H2 

e  +  H2 

elcistic 

0.0 

IffTlT-aa 

48 

e 

e  +  H 

— 

e  +  H 

elastic 

0.0 

4.7  X  10-^® 

1.6  X  10-1® 

36 

i 

H++H2 

— 

H+  +  H2 

elastic 

0.0 

7.9  X  10-1® 

49 

i 

H+  +  H 

— 

H+  +  H 

elastic 

0.0 

1.1  X  10-1^ 

Table  2.  The  Seven  Reactions  Used  in  the  Numerical  Model. 


Table  2  also  gives  representative  reaction  rates  for  the  inelastic  chemical  reactions  and 
representative  cross  sections  <rn  for  the  elastic  scattering  reactions.  These  representative  rates  and 
cross  sections  are  those  that  the  model  calculates  for  temperatures  of  leV  and  10 eU. 

Because  the  model  reduces  the  number  of  interacting  species  to  four  and  the  number  of 
interactions  to  seven,  most  of  the  complexity  of  the  hydrogen  ion  source  is  lost.  Unfortunately,  the 
dominant  reactions  in  a  hydrogen  plasma  at  the  energies  of  interest  are  among  those  discarded. 
The  model  does  not  consider  the  dominant  path  of  ionization  at  the  electron  temperature  that 
it  calculates,  nor  does  it  consider  the  dominant  paths  for  recombination.  As  a  result,  the  model 
calculates  a  charged  particle  balance  between  production  and  outflow  instead  of  the  correct  balance 
between  production  and  loss. 
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3.J  Ionization 


The  model  considers  only  the  ionization  reaction  between  electrons  and  neutral  hydrogen 
atoms.  Its  standard  initial  conditions  [9:422]  assume  that  the  density  of  neutral  molecules  is  larger 
than  the  density  of  neutral  atoms  by  over  an  order  of  magnitude.  Since  the  reaction  rate  for  the 
ionization  of  neutral  molecules  is  approximately  the  same  as  the  rate  for  the  ionization  of  neutral 
atoms  at  the  electron  temperature  calculated  by  the  model,  the  model  therefore  underestimates 
the  ionization  by  about  an  order  of  magnitude. 

Because  the  model  considers  chemical  reactions  among  only  the  species  H ,  H'*’,  and  e,  it 
produces  plasma  in  the  volume  by  way  of  Reaction  19  in  Table  2  and  in  Appendix  T: 

e  +  /f  — *  c  +  //■*■  +  e 

which  is  Reaction  2.1.5  in  Janev  ei  ai  [17:3].  Since  the  density  of  molecular  hydrogen  H2  is  roughly 
an  order  of  magnitude  larger  than  the  density  of  atomic  hydrogen  i/ ,  it  is  natural  to  consider  the 
contribution  of  the  ionization  reaction  between  electrons  and  molecular  hydrogen: 

e  +  H2  e  +  H2  +  e 

which  is  Reaction  24  in  Appendix  B  and  Reaction  2.2.9  in  Janev  ei  ai  [17:3].  Table  3  compares 
the  H  ionization  reaction  rate  ki^  for  Reaction  19  with  the  H2  ionization  reaction  rate  ^24  for 
Reaction  24  over  a  range  of  energies.  The  reaction  rates  given  for  Reaction  19  are  those  calculated 
by  the  model  [35],  which  agree  well  with  those  given  in  Janev  et  ai  [17:26],  and  the  reaction  rates 
for  Reaction  24  are  from  Janev  et  ai  [17:52]. 

Both  of  these  reaction  rates  show  a  sharp  increase  with  electron  temperature.  Because  the 
threshold  energy  for  Reaction  19  with  H  (cth  =  13.6  eV)  is  lower  than  the  threshold  energy  for 
Reaction  24  (cth  =  15.4  61^),  the  rate  of  Reaction  19  tends  to  dominate  somewhat  at  lower  tem- 
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kit  = 

)  (cm^/sec) 

T‘  (eV) 

ki9  {H) 

k24  (H2) 

1 

3.4  X  10-^^ 

<  10-^1 

3 

1.2  X  lO-io 

6.0  X  10"“ 

7 

2.7  X  10-® 

2.4  X  10-® 

10 

5.3  X  10-® 

6.3  X  10-® 

100 

3.2  X  10-^ 

5.0  X  10-^ 

Table  3. 


Comparison  of  Reaction  Rates  for  Ionization  Reaction  19  [35]  and  Ionization  Reaction  24 
[17:52]. 


peratures  (<  3  eV).  However,  at  the  electron  temperature  that  the  model  calculates  (~  7  eV), 
the  reaction  rates  for  both  reactions  are  approximately  equal,  so  the  higher  density  of  molecular 
hydrogen  causes  the  principal  ionization  path  to  be  Reaction  24.  By  this  measure,  the  model 
underestimates  the  rate  of  ionization  by  an  order  of  magnitude. 


S.2  Recombination 

Plasma  is  lost  in  the  volume  through  recombination  reactions  and  wall  losses.  The  primary 
electron  recombination  path  is  not  the  radiative  recombination  reaction  with  H*  (Reaction  26  in 
Table  2  and  Appendix  B)  upon  which  the  model  is  based,  but  rather  the  dissociative  recombina¬ 
tion  reactions  with  and  .  Because  these  reactions  have  multiple  products  to  carry  away 
momentum  and  energy,  they  are  fast  reactions.  On  the  other  hand.  Reaction  26  with  emits  a 
photon  to  conserve  energy  and  is  thus  a  slow  reaction  [34:18-19]. 

The  principal  electron  recombination  reaction  with  is  Reaction  2.2.14  in  Janev  [17:62]: 

where  H*  represents  an  electronically  excited  state  of  hydrogen.  Appendix  B  lists  a  less  specific  form 
of  this  dissociative  recombination  reaction  as  Reaction  16.  The  principal  recombination  reaction 
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with  has  two  possible  paths: 


H  +  H  +  H 

H2{v">h)  +  H'{n  =  2) 


which  are  Reactions  2.2.15a  and  2.2.15b  in  Janev  [17:64].  This  dissociative  recombination  reaction 
corresponds  roughly  to  Reaction  17  in  Appendix  B. 


Table  4  compares  the  reaction  rate  ibos  for  Reaction  26  with  the  reaction  rates  ki^  and  kn 
for  Reactions  16  and  17.  All  data  in  the  Table  4  are  estimated  from  the  curves  given  in  Janev  ei 
al.  [17].  The  table  shows  that  the  electron  recombination  rate  with  is  faster  from  1  to  3  eV , 


% 

II 

(cm^/sec) 

T‘{eV) 

kie  kn 

1 

1.4  X  10“^^ 

5.5  X  lO-**  4.2  X  10-“ 

3 

8.0  X  10-1^ 

2.8x10-®  2.6x10-® 

7 

4.7  X  lO-*'* 

1.7x10-®  2.0x10-® 

10 

3.6  X  lO"*'* 

1.3x10-®  1.5x10-® 

Table  4.  Reaction  Rates  [17]  for  Representative  Recombination  Reactions. 


and  the  rate  with  is  fastest  from  7  to  10  eV .  Either  path  for  neutralizing  electrons  by  way  of 
reactions  with  H}  or  is  more  probable  than  the  path  with  used  in  the  model  by  almost 
six  orders  of  magnitude.  Reaction  26  with  H'*'  is  clearly  inadequate  as  the  sole  loss  mechanism. 
The  result  of  the  lack  of  a  correct  loss  mechanism  is  that  the  model  achieves  an  incorrect  charged 
particle  balance  in  the  continuity  equation. 


3.3  The  Plasma  Balance 


The  plasma  balance  is  described  by  the  continuity  equation.  Equation  46  in  Appendix  A 
gives  a  form  of  the  continuity  equation: 


di 


~  ‘^production  ‘^loss 
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Since  the  final  solution  is  in  a  steady  state,  the  time  derivative  is  zero,  and  the  particle  balance 
along  the  axial  direction  becomes 


— (n"v;“)  =  n‘n°  ^{(Tnv)io„ization  -  n‘n'  recombinatio 


Under  the  assumptions  and  reduced  reaction  set  of  the  model,  this  balance  reduces  to 


dnVj 

dz 


nn°kig  —  (n)^k2e 


Using  the  model’s  published  standard  “test  case”  data  [9:422],  n  w  4  x  10“  cm“®,  n°  x  s» 
3  X  and  T*  as  7.5  eV.  At  T‘  =  7.5  eV,  the  model  computes  reaction  rates  of  kig  » 

3  X  10“®  cm“^scc“^  and  iS;26  as  5  x  10~^^cm~^scc“* ,  so  the  particle  balance  becomes 

^  ^  — ~ —  as  (4  X  10‘^)(3  X  10*^)(3  x  10“®)  -  (4  x  10^^)(4  x  10“)(5  x  10“^^) 

20  cm 

4  X  10^® cm“®sec"^  as  4  x  10*® cm“® sec“*  -8  x  10*°cm“^sec“* 

^  '  . . ^  ' - - - '  ' - 

flux  production  loss 

Thus  the  production  (ionization)  is  not  balanced  by  loss  (recombination)  as  would  be  ex¬ 
pected,  but  by  flow  out  of  the  extraction  region.  Since  all  of  the  plasma  that  the  source  produces 
is  extracted,  this  type  of  baletnce  would  be  ideal  for  an  ion  source;  however,  because  production 
is  underestimated  by  an  order  of  magnitude  and  because  loss  is  underestimated  by  several  orders 
of  magnitude,  it  is  not  realistic.  This  incorrect  particle  balance  is  a  significant  limitation  of  the 
model. 

In  simplifying  the  kinematics  of  the  hydrogen  ion  source,  the  model  discards  some  of  the 
dominant  reaction  paths.  It  underestimates  plasma  production  by  about  an  order  of  magnitude 
and  plasma  loss  by  almost  six  orders  of  magnitude.  The  resulting  incorrect  b  mce  between  pro¬ 
duction  and  outflow  instead  of  production  and  loss  significantly  limits  the  usefulness  of  the  model. 
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Section  4.4  shows  that  the  model  has  to  raise  the  electron  temperature  by  almost  an  order  of 
magnitude  to  make  up  for  the  underestimate  of  ionization. 
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IV.  The  Behavior  of  the  Unmodified  Model 


Running  the  model  with  the  published  inputs  for  the  test  rase  reproduced  the  results  in 
Glasser  and  Smith  [9:422]  as  accurately  as  could  be  determined  from  the  published  graphs.  However, 
the  calculated  behavior  of  the  primary  variables  differs  from  what  would  be  expected  from  intuition 
and  is  obtained  by  experiment.  Moreover,  results  differ  from  run  to  run,  even  when  the  initial  inputs 
are  identical. 

4-1  Replication  of  the  Results  in  the  Original  Article 

Figure  11  shows  the  temporal  evolution  of  the  four  primary  variables  (charged  particle  density, 
drift  velocity,  electron  temperature,  and  ion  temperature).  These  results  reproduce  those  in  Figures 
3  through  6  in  Glasser  and  Smith  [9:422].  Each  graph  displays  the  value  of  one  of  the  primary 
variables  as  a  function  of  the  axial  dimension  of  the  ion  source.  The  graphics  routine  plots  one 
curve  at  each  time  cycle  (—  3xl0~®  sec),  with  25  cycles  displayed.  Clustering  of  the  curves  indicates 
convergence  to  an  equilibrium  solution. 

Figure  12  reproduces  the  results  of  Figure  7  in  Glasser  and  Smith  [9:423],  showing  the  move¬ 
ment  of  the  node  positions  used  in  the  method  of  moving  finite  elements  as  functions  of  time.  The 
horizontal  axis  plots  time,  and  the  vertical  axis  plots  the  axial  dimension  of  the  source.  Node 
positions  tend  to  cluster  where  the  second  derivatives  of  the  primary  variables  are  the  greatest, 
which  occurs  in  this  case  at  inflection  points  on  either  side  of  the  magnetic  filter  field.  Convergence 
to  the  final  positions  is  rapid,  being  essentially  complete  by  the  seventh  cycle. 

The  time  order  of  the  plots  in  Figure  11  isn’t  completely  obvious.  For  this  reason,  Figure  13 
shows  the  values  of  the  primary  variables  at  the  outflow  boundary  (position  z  =  25  cm,  integration 
bin  j  =  40).  These  graphs  provide  a  means  to  verify  the  convergence  of  solutions  and  to  identify 
other  phenomena,  such  as  limit  cycles. 
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4-2  Plasma  Density 


The  upper  left  graph  in  Figure  11  shows  multiple  plots  of  the  calculated  values  of  the  plasma 
density,  n,  as  a  function  of  the  axial  dimension  of  the  ion  source.  The  lower  graph  in  Figure  14  shows 
the  equilibrium  state  of  the  plasma  density  after  25  integration  cycles  (about  8  x  10“'*  sec).  The 
plctsma  density  is  greatest  in  the  extraction  region.  This  result  seems  counterintuitive:  since  in  a  real 
ion  source  the  magnetic  filter  field  acts  as  a  barrier  to  ionizing  fast  electrons,  confining  them  to  the 
production  region,  and  since  electrons  and  ions  recombine  in  and  flow  out  of  the  extraction  region, 
one  would  expect  the  density  in  the  extraction  region  to  be  /otcer  than  in  the  production  region, 
rather  than  higher.  Moreover,  the  plasma  flows  along  the  axis  from  the  production  region  to  the 
extraction  region,  and  particles  should  flow  from  high  to  low  density,  rather  than  from  low  to  high. 
Experimentally,  according  to  Ehlers  and  Leung,  “the  plasma  density  in  the  extraction  chamber  is 
always  less  than  in  that  in  the  source  chamber”  [8:1423].  Johnson  [18:162]  gives  experimental  values 
of  the  axial  variation  of  plasma  density  that  confirm  the  observation  of  Ehlers  and  Leung.  The 
upper  graph  in  Figure  14  reproduces  Johnson’s  data  for  an  ion  source  with  a  magnetic  filter  field  of 
strength  80  Gauss  centered  at  15  cm.  Apparently  the  initial  guess  at  the  plasma  density  in  the  input 
to  the  test  case,  shown  in  the  graph  of  density  in  Figure  11,  is  roughly  correct;  however,  the  model 
takes  this  correct  solution  and  immediately  distorts  it,  as  early  as  the  first  cycle  (~  3  x  10"®  sec), 
into  one  very  similar  to  the  final  output. 

Another  troubling  result  is  that  the  final  values  of  the  plasma  density  depend  upon  the  value 
of  the  initial  guess  used  in  its  calculation.  This  is  troubling  because  there  is  no  physical  significance 
to  the  initial  guess  used.  The  path  to  the  final  answer  should  reflect  the  character  of  the  initial  guess, 
but  the  final  answer  should  depend  only  on  the  actual  physical  parameters.  After  this  behavior 
became  apparent,  the  model  was  run  with  several  different  combinations  of  initial  densities  at 
the  right  and  left  boundaries.  Three  initial  densities  were  chosen  to  match  the  “test  case”  initial 
conditions  in  the  original  article  [9:422]:  4.0  x  10**,  2.5  x  10**,  and  1.2  x  10**  cm~^  as  the  respective 
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Figure  14.  Axial  Variation  of  Plasma  Density:  Experimental  Measurement  (Top)  from  Johnson 
[18:162]  Compared  to  Calculations  by  the  Model  (Bottom).  In  both  graphs,  the  fil¬ 
ament  cathodes  are  located  at  r  =  0,  and  vertical  lines  mark  the  1/e  points  of  the 
magnetic  filter  field. 
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high,  medium,  and  low  densities.  The  high  density  was  the  left  initial  density  (at  the  filaments) 
in  the  standard  run,  the  low  density  was  the  right  initial  density  (at  extraction),  and  the  medium 
density  was  roughly  in  the  middle.  The  model  was  run  with  flat  initial  density  profiles  at  all  three 
values,  with  the  standard  initial  profile,  and  with  a  reversed  initial  profile.  Table  5  gives  the  results 
for  the  meiximum  final  density  resulting  from  these  runs  of  the  model,  grouped  by  the  initial  density 
at  the  left  boundary: 


initial  density 

initial  density 

maximum 

at  left  boundary 

at  right  boundary 

final  density 

4.0  X  10“ 

4.0  X  10“ 

5.326  X  10“ 

4.0  X  10“ 

1.2  X  10“ 

5.867  X  10“ 

2.5  X  10“ 

2.5  X  10“ 

3.466  X  10“ 

2.5  X  10“ 

1.2  X  10“ 

3.406  X  10“ 

1.2  X  10“ 

4.0  X  10“ 

1.645  X  10“ 

1.2  X  10“ 

1.2  X  10“ 

1.143  X  10“ 

Table  5.  Variation  of  the  Maximum  Final  Values  of  the  Plasma  Density  with  Changes  to  the 
Initial  Values. 


Table  5  shows  that  the  final  value  of  the  density  depends  strongly  upon  the  value  assumed 
for  the  density  at  the  left  boundary,  nj,  and  hardly  at  all  upon  the  value  at  the  right  boundary, 
n%max-  Section  5.1  shows  that  this  dependence  occurs  because  of  the  way  that  the  model  defines 
the  energy  flux. 

Another  feature  of  the  calculated  density  profile  is  the  apparently  parabolic  decrease  in  density- 
in  the  production  region.  It  is  possible  to  derive  this  parabolic  decrease  analytically.  Beginning 
with  the  one-dimensional  momentum  equation  implemented  in  the  model.  Equation  9  in  Section  6.1 
and  Equation  3.4  in  Glasser  and  Smith  [9:412] 


p®  dn 

4 

m  ^ 

a 

(n)2 

3m(n] 

-  -n“r“ 

[n  1' 

3*' 

az )] 

^  fix 
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and  discarding  all  but  the  two  dominant  terms,  Term  3  and  Term  9  [35]  identifies  an  approximate 
balance  between  these  two  terms: 


1  dn 

m  (n)2  dz 


+  dV,=0 


dfx 

Substituting  p“  =  nT“  and  solving  for  —  gives 


52  -  =  __!25_r 


where  F  =  nV^  is  the  particle  flux.  Figure  15  shows  that  the  model  calculates  a  constant  tem¬ 
perature  in  the  production  region,  and  Figure  19  shows  that  it  calculates  a  flux  in  the  production 
region  that  varies  approximately  linearly  with  position.  Denoting  the  slope  of  the  graph  of  the  flux 
by  A,  the  relation  r  =  Az  then  describes  this  linear  variation.  Substituting  this  linear  relation  for 
the  flux  into  the  simplified  momentum  equation  gives  a  differential  equation  in  z: 


dn  Amd 


Integrating  this  differential  equation  with  respect  to  z  yields  an  explicit  expression  for  the  plasma 
density: 


Amd 


where  zq  is  the  constant  of  integration.  Define  the  constant  cq  to  be  the  coefficient  of  the  (z)^  term; 
then  the  expression  for  the  density  takes  the  form 


n  =  zo-  co(?)^ 


which  describes  a  parabolic  decrease  in  the  plasma  density  as  the  distance  from  the  filament  cathode 
increases.  Figure  14  shows  a  similar  but  larger  decrease  in  the  plasma  density  in  the  experimental 
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values  from  Johnson  [18:162].  However,  Chapter  III  shows  that  the  plasma  balance  in  the  model 
does  not  correspond  to  that  in  a  real  source.  For  this  reason,  the  assumptions  of  the  derivation 
above  may  not  apply,  and  the  density  profile  in  the  production  region  of  a  real  source  may  decrease 
for  a  different  reason. 

Therefore,  even  though  the  model  calculates  a  density  profile  that  appears  to  be  qualitatively 
correct  in  the  production  region,  it  calculates  a  profile  through  the  magnetic  filter  field  that  is 
opposite  to  that  expected  from  intuition  and  observed  by  experiment.  Moreover,  the  final  values 
of  the  density  depend  on  the  initial  guess  at  a  solution.  Both  of  these  features  severely  limit  the 
utility  of  the  model  in  predicting  the  behavior  of  an  ion  source. 

4.3  Drift  Velocity 

The  upper  right  graph  in  Figure  11  shows  multiple  plots  of  the  time  evolution  of  calculated 
values  for  the  velocity  as  a  function  of  the  axial  dimension  of  the  ion  source.  Since  the  model 
assumes  that  the  plasma  is  initially  stationary,  the  initial  (zeroth-cycle)  plot  of  the  drift  velocity  in 
Figure  11  is  uniformly  zero  across  the  bottom  of  the  graph.  The  second  (first-cycle)  plot  of  the  drift 
velocity  overshoots  the  eventual  equilibrium  value,  to  which  it  converges  rapidly.  By  the  seventh 
cycle,  subsequent  curves  cannot  be  distinguished  from  one  another  on  the  scale  of  the  plot.  The 
upper  right  graph  of  Figure  13  shows  this  overshooting  and  rapid  convergence  for  the  drift  velocity 
at  the  outflow  boundary. 

The  average  drift  velocity  is  held  to  zero  at  ?  =  0  by  a  boundary  condition  applied  in  the 
module  arrays.  The  physical  interpretation  of  the  plot  is  as  follows:  Plasma  is  “born”  at  rest 
in  the  vicinity  of  the  filaments  at  z  =  0.  The  linear  increase  in  velocity  in  the  driver  region  in 
Figure  11  shows  that  the  acceleration  of  the  plasma  is  constant  there.  The  magnetic  filter  retards 
the  electrons  as  they  diffuse  across  it,  causing  a  reduction  in  the  average  drift  velocity.  After  the 
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filter  field  the  acceleration  once  again  becomes  constant.  This  behavior  appears  to  be  intuitively 
reasonable. 

4.4  Electron  Temperature 

The  lower  left  graph  in  Figure  11  shows  the  evolution  of  the  electron  temperature  T‘  as  a. 
function  of  the  ajcial  dimension  of  the  ion  source.  The  initial  (zeroth-cycle)  plot  of  the  electron 
temperature  is  near  the  bottom  of  the  graph,  and  the  second  (first-cycle)  plot  is  at  the  top.  The 
lower  left  graph  in  Figure  13  shows  that  the  second  calculated  value  of  the  electron  temperature  at 
the  outflow  boundary  overshoots  the  equilibrium  value  and  that  subsequent  calculations  converge 
to  it,  although  not  quite  as  rapidly  as  in  the  case  of  the  drift  velocity. 

The  lower  graph  of  Figure  15  shows  the  equilibrium  values  of  the  electron  temperature  at 
the  25‘^  cycle.  The  electron  temperature  decreases  across  the  magnetic  filter,  as  expected  from  the 
discussion  in  Section  2.3;  however,  there  is  a  problem  here,  too.  The  average  electron  temperature 
determined  experimentally  in  a  real  ion  source  is  approximately  1  eV  in  the  extraction  region 
[1:21]  and  2  eV'  in  the  production  region  [34:11],  not  6  eV  and  7  eV  as  predicted  by  the  model. 
Figure  15  compares  the  electron  temperature  predicted  by  the  model  on  the  lower  graph  with  the 
experimental  values  determined  by  Johnson  [18:162]  on  the  upper  graph. 

This  discrepancy  in  electron  temperature  can  be  explained  from  the  assumptions  in  the  model. 
The  model  assumes  one  single-temperature  maxwellian  distribution  of  electrons.  In  reality  the  dis¬ 
tribution  is  non-maxwellian,  but  it  can  be  described  reasonably  accurately  as  a  superposition  of 
two  maxwellian  distributions,  each  with  its  own  temperature.  As  described  in  Section  2.1,  the 
larger  of  the  two  contains  low-energy  “secondary”  electrons,  and  the  smaller  contains  high-energy 
“primary”  electrons.  Figure  16  shows  how  the  seemingly  small  discrepancy  of  ignoring  the  smaller 
distribution  can  make  an  order  of  magnitude  difference  in  the  results.  The  solid  lines  represent  the 
two-population  approximation  to  the  experimentally-determined  non-maxwellian  electron  distribu- 
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Figure  15.  Axial  Variation  of  (Secondary)  Electron  Temperature:  Experimental  Measurement 
(Top)  from  Johnson  [18:162]  Compared  to  Calculations  b\  the  Model  (Bottom).  In 
both  graphs,  the  filament  cathodes  are  located  at  z  =  0,  and  vertical  lines  mark  the 
1/e  points  of  the  magnetic  filter  field. 
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Figure  16.  Schematic  of  the  Difference  Between  a  One- Population  and  a  Two- Population  Model 
of  the  Electron  Energy  Distribution. 


tion.  The  dashed  line  represents  the  one-population  distribution  assumed  in  the  model.  Only  those 
electrons  with  energies  greater  than  the  threshold  energy  for  ionization  can  ionize  neutral  atoms  or 
molecules.  The  probability  of  an  ionization  is  proportional  to  the  ionization  cross  section  at  that 
energy.  Thus  high  energy  electrons  (the  shaded  area  in  the  diagram)  accomplish  most  of  the  ioniza¬ 
tion  in  the  tandem  volume  source.  By  using  a  single  electron  distribution,  the  model  requires  the 
electron  distribution  to  be  at  a  much  higher  temperature  to  achieve  a  given  rate  of  ionization  than 
it  would  if  it  were  to  use  two  distributions  at  different  temperatures.  This  is  a  serious  limitation  in 
the  model  which  should  be  corrected  if  the  model  is  to  correctly  predict  experimental  results. 


4.5  Jon  Temperature 

The  lower  right  graph  in  Figure  11  shows  multiple  plots  of  the  evolution  of  the  axial  depen¬ 
dence  of  the  ion  temperature  T’  .  Since  the  mode!  assumes  a  constant  ion  temperature  throughout 
the  source,  the  initial  (zeroth-cycle)  plot  of  the  ion  temperature  is  flat  across  the  middle  of  the 
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graph.  The  first-cycle  plot  transforms  this  flat  profile  into  a  form  very  much  like  the  equilibrium 
result.  The  lower  right  graph  in  Figure  13,  of  the  ion  temperature  at  ‘he  outflow  boundary,  shows 
that  the  ion  temperature  decreases  slightly  in  the  second  cycle  and  then  increases  until  it  overshoots 
and  then  settles  to  its  equilibrium  value.  The  equilibrium  profile  remains  relatively  flat,  varying 
only  by  about  six  percent  from  its  highest  to  its  lowest  value. 

In  contreist  to  its  prediction  for  the  electron  temperature,  the  model  predicts  an  ion  tempera¬ 
ture  T‘  that  is  slightly  higher  in  the  extraction  region  than  in  the  production  region.  Since  in  a  real 
volume  source,  ions  are  heated  by  collisions  with  hot  electrons  and  cooled  by  collisions  with  neutral 
gas  molecules  [18],  one  would  expect  the  opposite  behavior:  In  the  production  region,  the  ions  can 
collide  with  hotter  electrons  and  can  thus  receive  more  thermal  energy  there  than  in  the  extraction 
region.  Johnson  assumes  that  the  positive  ion  temperature  is  approximately  equal  to  the  negative 
ion  temperature  (T*  «  ),  which  he  calculates  to  be  within  the  range  0.2  61^  <  <161^ 

[18:263].  The  lower  right  graph  in  Figure  11  shows  that  the  ion  temperature  calculated  by  the 
model  falls  within  this  range. 

4-6  Results  of  Multiple  Runs  of  the  Code  (Irreproducible  Results) 

The  most  distressing  feature  that  the  model  exhibits  is  that  running  it  with  exactly  the  same 
initial  conditions  does  not  always  produce  the  same  results.  Table  6  shows  the  results  of  multiple 
runs  of  the  model  with  the  same  inputs  as  the  standard  case  published  in  Glasser  and  Smith  [9:422]. 
The  order  of  the  listing  runs  from  lowest  meiximum  density  to  highest. 

As  seen  in  Table  6,  the  runs  with  the  lowest  and  highest  maximum  density  are  greatly  different 
from  the  others.  Although  the  other  runs  have  results  similar  to  each  other,  there  is  still  a  significant 
variation  among  them. 

After  this  variation  from  run  to  run  was  observed,  the  standard  case  was  run  three  times  with 
enhanced  graphics.  Two  of  the  runs  were  similar,  yet  not  exactly  the  same.  The  plasma  densities 
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run 

^max 

^2  )max 

T* 

max 

1 
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X 
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X 

10^ 

X 

10^ 

X 

10- 

-1 

2 

5.104 

X 

SQ 

1.848 

X 

10® 
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X 

10° 
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X 

10- 

-1 

3 

5.120 

X 
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1.864 

X 

10® 

X 

10° 

X 

10- 

-1 

4 

5.138 

X 

B 

1.843 

X 

10® 

7.567 

X 

10° 

3.172 

X 

10- 

-1 

5 

5.149 

X 

9 

1.844 

X 

10® 

7.548 

X 

10° 

3.172 

X 

10- 

-1 

6 

5.163 

X 

i  9 

1.843 

X 

10® 

X 

10° 

3.173 

X 

10- 

-1 

7 

5.173 

X 

,9 

1.839 

X 

10® 

7.554 

X 

10° 

3.173 

X 

10- 

-1 

8 

5.840 

X 

HiBl 

1.615 

X 

10® 

7.276 

X 

10° 

X 

10- 

-1 

Table  6.  Irreproducible  Results:  Variation  of  Final  Maximum  Values  of  the  Primary  Variables 
among  Eight  Runs  of  the  Model. 

differed  in  the  fourth  decimal  place  after  the  first  cycle  (~  3  x  10~®  sec),  but  by  the  25*^  cycle 
(~  8  X  10“^  sec)  they  differed  in  the  second  decimal  place.  One  of  these  two  runs  is  given  as  the 
replica  of  the  test  case  in  Figures  11,  12,  and  13. 


The  third  run  does  not  even  converge  by  the  25‘**  cycle.  Figure  17  shows  the  results  of  this 
curious  run,  which  was  launched  from  the  identical  conditions  that  generated  Figures  11  and  13. 
This  figure  corresponds  to  Figure  13  in  that  it  shows  graphs  of  the  values  of  the  primary  variables 
at  the  outflow  boundary  as  functions  of  the  number  of  integration  cycles.  The  plasma  density 
appears  to  converge  normally  until  about  the  18***  cycle  (~  6  x  10““*  sec),  when  it  starts  to  diverge. 


To  verify  that  the  added  graphics  modules  did  not  cause  the  variation  from  run  to  run,  the 
model  was  run  three  times  without  them.  These  three  runs  also  exhibited  similar  variations  from 
run  to  run. 


This  variation  from  run  to  run  is  the  most  serious  limitation  of  the  model.  If  the  results  are 
not  reproducible,  then  the  model  is  practically  useless.  Any  further  research  into  the  behavior  of 
the  model  should  idenvtify  and  correct  this  limitation. 


4-7  Summary  of  the  Results  of  the  Unmodified  Model 


In  summary,  even  though  the  model  was  able  to  approximately  reproduce  the  published 


results  in  Glasser  and  Smith  [9:422],  it  is  limited  by  several  factors.  First,  calculated  values  of  the 
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plasma  density  trend  opposite  to  intuition  and  experiment.  Second,  the  final  values  of  the  density 
depend  on  the  initial  guess  in  the  input  file.  Third,  the  calculated  electron  temperature  is  almost 
an  order  of  magnitude  higher  than  experimentally  observed.  Finally,  and  most  seriously,  calculated 
results  differ  run  to  run,  even  when  the  initial  inputs  are  identical. 
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V.  The  Boundary  Energy  Flux 


This  thesis  takes  two  exceptions  with  the  approach  that  the  model  uses  to  establish  and 
distribute  the  energy  input  from  the  discharge.  First,  the  model  determines  the  energy  flux  at  the 
boundary  of  the  driver  region  from  the  initial  guess  at  the  correct  values  of  the  density  and  the 
temperature.  Thus  the  final  solution  also  depends  on  these  initial  conditions,  which  are  not  physical 
parameters  of  the  system.  Second,  the  model  injects  energy  from  the  discharge  into  only  the  first 
integration  bin  at  z  =  0.  This  leftmost  integration  bin  is  where  the  model  locates  the  cathode 
filaments  (see  Figure  3  in  Chapter  I).  Because  the  mean  free  path  of  the  primary  electrons  is  of 
the  order  of  the  length  of  the  source,  a  uniform  deposition  of  discharge  energy  into  the  production 
region,  with  the  energy  flux  at  the  boundary  set  approximately  to  zero,  would  appear  to  be  a  better 
representation. 

5.J  Dependence  of  Density  on  the  Left  Initial  Condition 

Table  5  in  Section  4.2  shows  that  the  final  values  of  the  plasma  density  depend  strongly  upon 
the  initial  value  assumed  for  the  density  at  the  left  boundary  (z  =  0).  This  dependence  occurs 
because  of  the  way  that  the  model  defines  the  energy  flux  at  the  left  boundary.  In  line  174  of  the 
main  module  of  the  pos  code,  the  definition  for  the  left  boundary  energy  flux  efiux  is; 

174  •flux=-«llfac*chiprpe(0)*(yr(3)-yl(3))/xinax 

This  left  boundary  energy  flux  represents  the  flow  of  energy  density  from  the  electrical  dis¬ 
charge  into  the  plasma.  In  symbols,  the  expression  for  the  left  boundary  energy  flux  F'-g  is; 

=  ‘^0^z=0  ~ 

*mor 
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'j^e  _ 

where  ■  - iz_  jg  the  initial  electron  temperature  gradient,  is  the  initial  electron 

^max 

thermal  conductivity*  (at  the  left  boundary)  in  the  direction  perpendicular  to  the  magnetic  filter 
field,  and  Fq  is  a  “user  parameter”  that  specifies  how  much  energy  the  electrons  conduct  into  the 
plasma  [9:423].  The  input  file  to  the  test  case  sets  eflfac  (Fq)  to  0.02  [9:422]. 

Lines  117  and  118  of  the  module  brag  define  the  perpendicular  electron  thermal  conductivity 
chiprpe  as 

m 

where  p'  =  n*T®  defines  the  electron  pressure  p‘  ,  t‘  =  —  defines  the  electron  time  between 
collisions  r*  in  terms  of  the  electron  collision  frequency  i/‘  ,  x‘  =  defines  the  argument 
of  the  polynomial  V{x)  in  terms  of  the  electron  gyrofrequency  w*  and  the  mass  of  an  electron  m*. 
Thus  chiprpe(O)  is  or  the  electron  thermal  conductivity  /£"*■'  evaluated  at  the  left  boundary 

z  =  0. 

Substituting  these  values  into  the  expression  for  the  energy  flux  input  at  the  left  boundary 
and  rearranging  factors  shows  a  dependence  on  plasma  parameters  calculated  at  the  left  boundary 
at  <  =  0  and  on  initial  conditions  input  into  the  model: 


r;=.  =  - [ u.))  {".'..ir,. {Tt.,...  - TU,)) 

t  m  *mor  )  *>'»<<  ✓ 


constants 


plasma  parameters 


initial  conditions 


This  expression  shows  explicitly  that  the  energy  flux  input  at  the  left  boundary  depends  on  three 
initial  conditions. 


*  Braginskii  defines  the  thermal  conductivity  for  species  a  as 

«“  =  -i_(P(x“)) 


in  his  equations  4.37  luid  4.40  [3:249-250].  Glasser  and  Smith  set  k”  =  in  their  equations  3.17  and  3.18 

[9:413],  and  in  the  computer  code  they  use  chi  (x)  to  represent  the  thermal  conductivity  as  defined  by  Braginskii. 
This  thesis  follows  Braginskii 's  convention. 
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First,  the  energy  flux  input  depends  on  n*_o,  the  initial  density  at  the  left  boundary: 


r‘_o  oc  n*_o 

Table  5  in  Section  4.2  showed  that  the  final  solutions  the  model  calculates  depenr  on  the  initial  val¬ 
ues  of  the  density  at  the  left  boundary.  This  dependence  apparently  comes  through  the  calculation 
of  the  energy  flux  input  at  the  left  boundary.  Section  5.2  elaborates  on  this  observation. 

Second,  the  energy  flux  at  the  left  boundary  depends  on  T/_o,  the  initial  electron  temperature 
at  the  left  boundary: 

n=o «  Tu 

Third,  the  energy  flux  at  the  left  boundary  depends  on  difference  between  the  initial  electron 
temperature  at  the  left  boundary  and  the  initial  electron  temperature  at  the  right  boundary: 

The  last  two  dependences  on  the  initial  temperature  are  not  as  straightforward  to  identify  in  the 
output  from  the  model  as  is  the  first  dependence  on  the  initial  density  because  they  involve  not 
only  the  temperature  at  the  left  boundary  but  also  the  difference  in  temperature  between  the  two 
boundaries. 

5.2  Left  Boundary  Condition  for  the  Temperature 

The  only  place  that  the  computer  code  uses  the  variable  eflux  is  in  line  269  of  the  module 
arrays  (line  549  of  the  composite  code)  where  it  sets  the  left  boundary  temperature  gradient  yl  (3,  0) 
equal  to  the  ratio  of  the  energy  flux  and  the  thermal  conductivity  perpendicular  to  the  magnetic 
field: 
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269 


yl(3,0)=-eflux/chiprpe(0) 


or  in  symbols, 


dT‘\  n=o  ^ 

)  2  =  0  ^2  =  0 


_p£  ^2=0  ^£=ro)f-o 

°  2  «-■-«  \ 

Zmor  ^2=o)t=t> 


where  t'  is  the  current  time  of  integration. 


One  significance  of  this  boundary  condition  is  that  all  of  the  energy  injected  into  the  volume 
enters  the  first  bin  at  2  =  0.  Since  the  mean  free  path  of  an  electron  is  determined  primarily  by 
the  density  of  the  most  abundant  constituent,  H^,  then 


(4  X  10i3cm“^)(l  x  10“*®cm2)  4  x  10~2cm~'  25cm 

which  is  the  length  of  the  device.  Thus  first-generation  primary  electrons  can  interact  with  neutral 
gas  molecules  throughout  the  entire  volume  in  the  production  region,  so  the  model  should  input 
energy  from  the  discharge  uniformly  into  the  whole  production  volume,  not  just  into  the  first 
integration  bin  at  2  =  0. 


5.3  Summary 

This  chapter  h^ls  identified  two  reservations  about  the  method  that  the  model  uses  to  incor¬ 
porate  energy  input  from  the  electrical  discharge.  Because  the  model  bases  its  evaluation  of  the 
energy  flux  at  the  left  boundary  on  an  arbitrary  input  initial  guess  at  the  density  and  temperature 
at  the  boundaries,  and  because  it  maintains  this  value  of  the  flux  throughout  its  calculation,  the 
final  solution  also  depends  on  these  initial  conditions.  Moreover,  the  distribution  of  energy  input 
from  the  discharge  appears  to  be  inconsistent  with  the  electron  mean  free  path.  Because  the  mean 
free  path  of  the  priir~*y  electrons  is  of  the  order  of  the  length  of  the  source,  a  uniform  deposition 
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of  energy  into  the  production  region,  with  the  energy  flux  at  the  boundary  set  approximately  to 
zero,  appears  to  be  a  better  representation. 
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VI.  The  Derivation  of  the  Equations  of  the  Model 


The  basis  for  the  model  is  Braginskii’s  derivation  [3]  of  transport  equations  for  two  fully 
ionized  fluids  (designated  by  superscripts  i  and  e,  where  i  represents  H"^).  Glasser  and  Smith 
have  enhanced  these  transport  equations  with  additional  terms  to  simulate  the  effects  of  the  set  of 
chemical  reactions  listed  in  Table  2  in  Chapter  III  and  the  effects  of  collisions  with  neutral  species 
(H2  and  H)  [9:410],  The  pos  code  uses  these  equations  to  set  up  differential  equations  for  the  four 
principal  variables  (plasma  density,  drift  velocity,  electron  temperature,  and  ion  temperature)  to 
be  solved  by  the  dirk2  differential  equation  solver  [5], 

This  chapter  contains  a  general  overview  of  the  derivation  of  the  equations  of  the  model. 
Appendix  A  gives  a  detailed  derivation  of  the  three  velocity  moments  of  the  Boltzmann  Equation 
and  explains  the  symbols  used  in  greater  detail. 

6.t  Development  of  the  Transport  Equations 

Braginskii’s  derivation  of  the  transport  equations  for  two  fully-ionized  fluids  begins  with  the 
Boltzmann  Transport  Equation,  Equation  1.1  in  Braginskii  [3:205],  for  a  distribution  /“(v,x,t)  of 
a  particle  of  species  q: 


dr\ 


/  coll 


(4) 


where  n“  is  the  average  number  per  volume  of  particles  of  species  a;  x,  ,  v,  ,  and  a,  are  the  com- 
ponents  of  position,  velocity,  and  acceleration,  respectively;  and  1  is  the  effect  of  collisions. 

/coll 

Braginskii  successively  multiplies  Equation  4  by  1,  m°v°,  and  -m°v°v°  and  integrates  over  the 
velocity  to  obtain,  respectively,  the  particle  transport  equation,  the  momentum  transport  equation, 
and  the  energy  transport  equation  for  particles  of  species  a. 
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6.2  The  Continuity  Equation 


Integrating  Equation  4  over  the  velocity  gives  the  particle  transport  equation,  or  continuity 
equation,  Equation  1.11  in  Braginskii  [3:208],  which  is: 


dn°‘  _ 

dt  dxf 


(5) 


where  is  the  component  of  the  average  velocity  of  particles  of  species  a  in  the  direction 
and  Sfgii  is  the  net  production  of  charged  particles  (ionization  minus  recombination).  Thus  the 
continuity  equation  states  that  the  time  rate  of  change  of  the  number  density  of  charged  particles 
in  a  volume  of  space  is  equal  to  the  net  flow  of  particles  into  that  volume  plus  the  net  production 
of  particles  within  the  volume.  In  his  derivation,  Braginskii  assumes  no  production  or  annihilation 
of  charged  particles  in  the  fully-ionized  fluids  and  thus  sets  =  0.  On  the  other  hand,  the 
model  assumes  only  charge  quasineutrality,  n*  =  n'  =  n,  so  electrons  and  ions  are  produced  and 
annihilated  in  pairs  and  5*^,,,  =  ■^oll  =  Scou.  The  model  considers  production  (ionization)  from 
Reaction  19  and  loss  (recombination)  from  Reaction  26  in  Table  2  or  Appendix  B,  so  the  net 
production  is: 

Scoii  =  n*n°^•l9  -  n*n‘Ar26  =  rin°ki9  -  {n)^k26  (6) 

which  is  Equation  3.6  in  Glasser  and  Smith  [9:412]. 

The  model  assumes  ambipolar  drift,  Vf  =  V*  =  V^,  where  the  z-direction  is  along  the  axis  of 
the  ion  source,  with  z  =  0  at  the  filament  (cathode)  end  and  z  =  Zmax  at  the  extraction  end.  The 
one-dimensional  form  of  the  continuity  equation  is  thus: 


dn  d{nVi) 
dt  dz 


^coll 


(■) 
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6.3  The  Momentum  Equation 


Multiplying  Equation  4  by  m°‘vf  and  integrating  over  the  velocity  yields  the  momentum 
transport  equation,  or  momentum  balance  equation.  Equation  1.14  in  Braginskii  [3:209],  which  is: 


(8) 


where 

dv“ 

J  /  coll 

describes  the  momentum  exchange  through  collisions. 

From  Equation  8,  Glasser  and  Smith  [9:411]  derive  a  similar  one-dimensional  form  of  the 
momentum  transport  equation  as  their  Equation  3.2: 


dV  dV 

nm<^ +  ZenE,  ±  Rf  -  m“ V, Scoii  -  d" (9) 


Several  features  of  Equation  9  are  evident.  First,  Braginskii’s  has  been  taken  in  the  z 
direction  only  and  extended  to  include  the  effects  of  collisions  with  neutral  particles: 


=  ±Rf  4-  R^  - 


where  describes  the  exchange  of  momentum  in  the  z  direction  due  to  collisions  between 

charged  particles  of  species  i  and  e  (when  q  =  e  in  the  electron  equation,  then  0  =  j;  conversely, 
when  or  =  j  in  the  ion  equation,  then  0  =  e),  describes  the  exchange  of  momentum 

due  to  the  net  production  of  charged  particles  of  species  a  from  collisions  with  neutral  particles, 
d"l4  describes  the  exchange  of  momentum  due  to  collisions  of  particles  of  species  a  with  neutral 
particles  (the  drag  due  to  collisions  with  neutral  particles),  and  fZf"  appears  to  be  a  redundant 
momentum  exchange  term. 
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A  second  feature  of  Equation  9  is  that  it  is  missing  the  term  in  Equation  8  from  the  interaction 
of  the  velocity  with  the  magnetic  filter  field,  Let  the  z  direction  be  defined  along  the 

central  axis  of  the  volume  source,  from  the  driver  to  extraction  regions,  and  let  the  x  direction  be 
aligned  with  the  magnetic  filter  field  in  Figure  3  so  that  the  y  direction  is  perpendicular  to  both 
the  filter  field  and  the  2L\is.  Figure  18  depicts  this  geometry.  Since  the  filter  field  is  applied  in  one 
direction  only.  By  =  =  0.  Then  the  only  surviving  term  in  the  z  direction  of  the  interaction 

between  the  drift  velocity  and  the  filter  field  is  the  term  —VyB^.  The  model  assumes 

that  for  the  one-dimensional  case  Vj,  =  0.  Therefore,  the  only  effect  of  the  magnetic  filter  field 
enters  through  the  parameter  x“  defined  by 


where  w“  is  the  gyrofrequency  of  species  a  and  r“  is  the  mean  time  between  collisions.  The 
parameter  x“  appears  in  two  of  the  transport  coefficients:  first,  in  the  order  coefficients  of 
transport  of  viscosity,  77“  : 
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where  V(x°')  is  a  dimensionless  polynomial  in  x°‘,  and  second,  in  the  coefficients  of  thermal  transport 
in  the  direction,  /cf ; 


m“ 


The  one-dimensional  expression  in  Equation  9  should  correspond  to  similar  one-dimensional 
expressions  from  other  sources.  For  example,  Golant  [10:309]  gives  an  expression  for  steady-state 
electron  and  ion  motion  perpendicular  to  an  applied  magnetic  field  in  a  fully-ionized  plasma: 


T  eEx  -  ivj.  (nT")  q:  ^  [VJ  x  H]  q:  ^  [h  x  VT‘]  T  (V1  -  V^)  =  0  (10) 

where  the  upper  signs  are  for  the  electron  equation  (a  =  e)  and  the  lower  signs  for  the  ion  equation 
(a  =  i),  the  symbol  ±  represents  the  directions  perpendicular  to  the  magnetic  field,  the  magnetic 
field  intensity  H  is  related  to  the  magnetic  induction  B  by  pH  =  B,  with  p  representing  the 

magnetic  permeability,  P**  is  the  “frequency  of  electron-ion  collisions  averaged  over  the  Maxwellian 

Z^tH  H  B 

distribution”  [10:308],  is  the  cyclotron  frequency  for  species  o  [10:313],  and  h  =  —  =  — 

is  a  unit  vector  in  the  direction  of  the  magnetic  field  [10:315]. 

Under  the  condition  of  ambipolar  diffusion,  the  mean  velocities  of  the  electrons  and  the  ions 
are  approximately  equal  (V^  w  V]l),  so  the  last  term  of  Equation  10  is  zero.  If  the  x  eixis  lies  in 
the  direction  of  the  magnetic  field  as  in  Figure  18  so  that  By  =  Bt  =  0,  then  Golant’s  expression, 
in  indicial  notation,  becomes  two  analogous  equations,  one  for  motion  in  the  y  direction,  the  other 
for  motion  in  the  z  direction: 


^  ^  foT"!  I  \V°n  1  I  ^ 


1  d 


Ze 


3  ZpP' 


cp 


B 


-Br 


B, 


dT’ 

dz 

dy 


■‘y 


0 

0 
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where  Z  is  the  charge  state  (Z  =  1  for  ions  Z  =  — 1  for  electrons).  If  the  temperature  distribution  in 

ar* 

the  y  direction  is  uniform,  then  — —  =  0.  Upon  multiplying  by  n,  rearranging  terms,  and  applying 

oy 

the  ideal  gas  law  =  nT^ ,  the  equation  of  motion  in  the  z  direction  becomes: 


0  =  +  ZenEz  + 

az 


Ze 

cp 


(11) 


A  few  simplifying  assumptions  can  reduce  Equation  9  to  a  similar  form.  In  a  steady-state, 
the  first  term  is  zero.  Neglecting  the  smaller  terms  of  those  remaining.  Equation  9  becomes 


-  ■^p'^  +  ZenEz 

which  is  identical  to  Golant’s  form  in  Equation  11  except  for  the  missing  term  suggested 

by  Equation  8. 

In  reality,  the  magnetic  filter  doesn’t  go  away.  Electrons  orbit  in  the  filter,  and  ions  are 
deflected.  Both  processes  result  in  velocities  in  the  y  direction.  Thus  Vy  =  0  is  not  a  good 
assumption,  even  in  a  one-dimensional  model. 

The  third  feature  of  Equation  9  is  that  it  is  really  two  equations,  one  with  q  =  e  and  0  =  i, 
the  other  with  a  =  i  and  0  —  e.  Adding  these  two  equations  eliminates  the  electric  field,  Ez, 
and  the  coulomb  interaction  term,  The  resulting  equation  is  Equation  3.4  in  Glasser  and 

Smith  [9:412],  which  the  model  uses  as  the  basis  of  its  differential  equation  for  the  drift  velocity  Vz  '. 


dVz 

dt  ^  dz  m  ■"  (n)^  dz 


4  dVz 
3m(n)^  dz 


dn 

dz 


-I- 


-I-  —  i-RV'  -  RT 
nm 


+  rnVzScoii)  +  dVz  =  0 


(12) 


where  m  =  m'  m* ,  tj®  is  the  coefficient  of  viscosity  of  species  o,  r“  is  the  mean  time  between 
collisions  for  species  a,  and  of  is  a  composite  drag  defined  in  Equation  18. 
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6-4  The  Energy  Equation 


Multiplying  Equation  4  by  and  integrating  over  the  velocity  gives  the  energy  trans¬ 

port  equation,  which  is  Equation  1.23  in  Braginskii  [3:210]: 


i  (in-m-Vi-V;.  +  In-T-)  +  V-  + 

+  (n?.  V/)  +  qf}  =  -Zen'^EiV^  + 


(13) 


To  develop  differential  equations  for  the  temperatures  T*  and  T*  the  model  uses  a  one¬ 
dimensional  version  of  Equation  13  of  the  form: 


3  ar® 

5"-5r+ 


2  ^  dz  dz  m®  dz  \  dz  ) 


-  g“ +  im“(r)2s„«  (14) 


where  /c“  is  the  thermal  conductivity  of  species  a,  Q",„t  is  an  energy  “sink”  [9:412]  due  to  radiation 
from  and  chemical  reactions  involving  species  a,  and  g®  is  the  rate  of  heat  transfer  to  species  o, 
defined  as 

g®  =  g®'’  -I-  g®" 

where  g®^  is  the  rate  of  transfer  from  charged  species  0  to  charged  species  q,  and  g®"  is  the  rate 
of  transfer  from  neutral  particles  to  species  a. 


6.5  Differential  Equations 

The  computer  code  stores  each  of  the  primary  variables  n,  V^,  T',  T*,  and  z  £is  a  subarray  of 
j  elements  in  the  i**’  row  of  the  array  j/,y  (the  variable  y  (i,  j)  in  the  computer  code),  where  j  is  the 
number  of  the  integration  bin.  Thus  yij  =  nj,  j/2>  =  (f/)i,  yzj  -  Tj,  y^j  =  Tj,  and  ysj  =  Zj.  The 
model  solves  for  the  component  variables  and  isolates  the  time  dependence  of  the  one-dimensional 
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moment  equations  above  to  develop,  for  each  of  the  four  primary  variables  j/,  ,  differential  equations 
of  the  form  given  in  Equation  4.1  in  Glasser  and  Smith  [9:413]: 


dyi  _  c 
dt  ~  '  dz 

where  Si  is  a  matrix  of  source  terms,  and  Fi  is  a  matrix  of  flux  terms  defined  by  Equation  4.2  in 
Glasser  and  Smith  [9:413]: 

Fi  =  Ci  -  Di^ 

The  pos  code  then  exports  these  differential  equations  to  the  dirk2  code  that  solves  them. 


6.6  Numerical  Approach 

The  dtrkS  code  uses  a  diagonally  implicit  Runge-Kutta  method  of  order  2  and  the  method 
of  moving  finite  elements  to  solve  the  differential  equations  supplied  by  the  model.  The  method  of 
moving  finite  elements  concentrates  integration  “bins”  at  inflection  points  where  the  second  deriva¬ 
tives  of  the  primary  variables  are  the  greatest.  Thus  calculations  in  regions  where  the  gradients  of 
the  primary  variables  are  changing  sharply  use  a  finer  integration  mesh,  and  calculations  in  regions 
where  the  gradients  of  the  primary  variables  change  relatively  little  use  a  coarser  mesh.  By  means 
of  this  method,  the  model  achieves  greater  computational  accuracy  for  any  given  number  of  bins. 


6.7  Conclusion  to  Derivation  of  the  Equations 

Thus  Glasser  and  Smith  implement  Braginskii's  derivation,  from  the  velocity  moments  of  the 
Boltzmann  equation,  of  transport  equations  for  two  fully  ionized  fluids  in  one  dimension  as  the 
basis  for  the  model.  In  the  course  of  review  of  the  derivations,  three  inconsistencies  were  identified, 
where  the  equations  implemented  in  the  model  do  not  match  the  derivations.  These  inconsistencies 
are  the  subject  of  Chapter  VII. 
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VII.  Inconsistencies  in  the  Derivation  of  the  Equations 


There  are  three  inconsistencies  in  the  model  where  the  derivation  and  coding  of  the  differential 
equations  for  the  four  primary  variables  in  the  model  differ  from  the  one-dimensional  velocity 
moment  equations  adapted  from  Braginskii.  This  chapter  identifies  these  three  inconsistencies, 
describes  the  attempts  to  correct  them,  and  presents  the  results  of  those  attempts. 


7.1  Omitted  Momentum  Term  tn  the  Model 

The  first  inconsistency  is  that  in  lines  114,  127,  and  128  of  the  module  fluxes  the  computer 
code  omits  the  term  for  momentum  exchange  between  ions  and  neutral  particles,  from  the 
second  dissipative  source  term,  S^,  included  in  Equation  4.12  in  Glasser  and  Smith  [9:414]: 


^2 

n  az  m'n 


^  -m'vS.oii 

<  omitted  t 


(15) 


7.2  Extra  Temperature  Gradient  in  the  Model 

The  second  inconsistency  is  that  in  setting  up  the  differential  equations  for  temperature, 
the  implementation  of  the  computer  code  squares  the  temperature  gradient  in  one  term  of  the 
energy  transport  equation.  Instead  of  the  correct  form  given  by  Equation  14,  the  equation  actually 
implemented  in  the  computer  code  is: 
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(16) 


Dimensional  analysis  shows  that  the  form  of  the  fourth  term  is  correct  in  Equation  14.  not 
in  Equation  16.  If  [f]  represents  a  dimension  of  length,  [m]  of  mass,  and  [t]  of  time,  then  the  first 
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term  of  Equations  14  and  16  shows  that  they  have  dimensions  of  power  density: 


[3 


2 


dt 


[mm?] 

[t] 


The  fourth  term  in  Equation  14  also  has  dimensions  of  power  density; 


'  1  d 
m"  dz 


[^]w[^]  [4]  >'> 


meyt- 


as  required.  Thus  multiplication  by  the  “added”  gradient  of  the  temperature  in  Equation  16  is 
incorrect. 


7.3  Extra  Velocity  Gradient  in  the  Model 

The  third  inconsistency  occurs  in  the  derivation  of  the  differential  equation  for  14  from  the  sum 
of  the  one-dimensional  momentum  equations.  Instead  of  implementing  the  momentum  equation 
in  the  form  of  Equation  12,  the  computer  code  implements  it  with  the  gradient  of  the  velocity 
squared: 
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Once  again  dimensional  analysis  shows  that  the  form  developed  in  Section  6.1  as  Equation  12 
is  correct.  The  dimensions  of  the  first  terms  in  Equations  12  and  17  are  dimensions  of  acceleration: 
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The  leading  terms  in  front  of  the  parentheses  in  the  fifth  term  in  Equations  12  and  17  have  dimen¬ 
sions  of 


-l-E- 

'  1  ■ 

m  dz  ^  n 
.  01 

m 

m  = 


which  are  also  dimensions  of  acceleration,  so  the  quantity  inside  parentheses  in  the  fifth  term 

must  be  dimensionless.  There  is  also  a  simpler  way  to  see  this.  For  the  two  expressions  inside 

the  parentheses  to  be  subtracted,  they  must  have  the  same  dimensions.  Since  1  is  dimensionless, 
4 

must  be  dimensionless,  too.  Its  dimensions  are 

3  dz 
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Thus  the  correct  dimensionless  term  is  the  one  without  the  extra  gradient  of  velocity. 


7.4  Attempts  to  Correct  Inconsistencies  in  the  Equations 


Correcting  the  first  inconsistency  by  adding  the  missing  f?'"  term  (sc  (3,  j)  in  the  computer 
code)  back  into  the  equation  for  (dsource  (2,  j)  in  the  computer  code)  is  simple.  The  modification 
changed  lines  127  and  128  of  the  module  fluxes  from 


127  dsourc«(2, j)-d80urce(2, j) 

128  2  +(sc(2, j)-mi*yy(2, j)*sc(l. j))/(mi'»yy(l, j)) 


to 


127  dsourc«(2, j)=d80urce(2,j) 

128  2  +(8c(2, j)+8c(3, j)-ini*yy(2, j)*8c(l, j))/(mi*yy(l , j)) 


where  the  sc  (3,  j)  term  has  been  added  to  the  continuation  line. 


Correcting  the  second  and  third  inconsistencies  was  a  more  involved  matter  because  the  code 
divides  the  equations  up  differently  from  the  way  the  derivations  do.  Simply  dividing  through  by 
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the  extra  terms  where  they  appear  resulted  in  “divide  by  zero”  errors.  For  example,  the  first  call  of 
the  subroutine  fluxes  occurs  while  the  drift  velocity  is  still  initialized  to  zero,  so  the  velocity  cannot 
be  a  divisor.  This  problem  was  avoided  by  not  dividing  by  the  terms  to  be  removed;  instead,  the 
extra  terms  were  removed  from  variables  where  they  were  multiplicative  factors.  So  eis  not  to  alter 
other  expressions  derived  from  these  variables,  the  terms  had  to  multiplied  back  in  at  other  places. 

For  example,  to  correct  the  second  inconsistency,  removing  the  extra  temperature  gradients 
from  the  energy  equations  required  changing  lines  109  through  110  and  115  through  118  of  the 
fluxes  module.  The  unmodified  code  read: 


109  qte=-chiprpa(j)*boltz*tel 

110  qti=-chiprpi(  j)’»boltz*til 

116  dsourceO,  j)=-(boltz*(qeilac(j)*(te-ti)+qeOfac(j)*(te-teinpO)) 

116  2  +pizze*vl+qtea(denl/den))/(1.6*boltz*den) 

117  dsource(4, j)=-(boltz*(qeilac( j)*(ti-te)+qi01ac(j)*(ti-tenip0)) 

118  2  +pizzi*vl+qti*(denl/den))/{1.6*boltz*den) 


The  modification  removes  the  two  gradient  terms  tel  and  til  from  qie  and  qti  in  lines  109 
and  110  and  adds  them  back  into  dsource  (3,  j)  and  dsource  (4,  j)  in  lines  115  through  118. 


109  qte=-chiprpo(j)*boltz 

110  qti=-chiprpi(j )*boltz 

115  dsourceO, j)=-(boltz*(qeilac( j)*(ta-ti)+qe01ac(j)*(te-temp0)) 

116  2  +pizze*vl+(qte*tel)*(denl/den))/(1.6*boltz*den) 

117  dsource(4, j)=-(boltz*(qeilac(j)*(ti-te)+qiCiac(j)*(ti-tempO) ) 

118  2  +pizzi*vl+(qti*til)*(denl/den))/(l .5*boltz*den) 


The  modification  had  the  same  effect  as  dividing  the  terms  dflux  (3,  j)  and  dflux  (4,  j)  in 
lines  112  and  113  by  tel  and  ieS  respectively,  but  without  the  “divide  by  zero”  error  in  the  initial 
call  of  fluxes. 

.  9v 

The  correction  of  the  third  inconsistency  required  dividing  the  gradient  of  the  velocity  — 
(I’l  in  the  computer  code)  out  of  Z?2  {dflux  (2,  j)  in  the  computer  code).  To  avoid  the  error  on 
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the  initial  call,  vl  was  removed  from  pizze  (11'^)  and  pizzi  (11^^)  in  lines  106  and  107  of  fluxes.  To 
compensate  for  removing  vl  from  the  Ilf^  terms,  vl  was  multiplied  back  into  dsource  (2,  j)  (5^), 
dsource  (3,  j)  (S3),  and  dsource  (4,  j)  (S4).  The  original  code  was: 


106  pizzea-fourby3*(etae0(j)/3.+etael(j))*vl 

107  pizzis-f  ourby3*  (etaiO  (  j  )  /3 .  -t-atail  (  j  )  )  *vl 

108  pizz=pizze-tpizzi 

109  qte=-chiprpe( j)*boltz*tel 

110  qti=-chiprpi( j)*boltz*til 

111  dilux(2, j)=pizz/(den*mi) 

112  dllux(3, j)=qt«/(l .6*boltz*den) 

113  dllux(4, j)=qti/(l .6*boltz*den) 

114  dsource ( 2 , j ) =-df lux ( 2 , j ) * (deni /den) -drag ( j ) *v 

116  d80urce(3, j)=-(boltze(qeiiac(j)e(te-ti)+qeOfac(j)*{te-tempO)) 

116  2  +pizzeevl+qte*(denl/den))/(1.5*boltz*den) 

117  dsource(4, j)=-(boltz*(qeilac(j)e(ti-te)+qi01ac(j)*(ti-temp0)) 

118  2  +pizzi’»vl+qti*(denl/den))/(1.5*boltz*den) 


After  removing  vl  from  dflux  (2,  j),  the  modified  code  was: 


106  pizzes-f  ourby3e  (  etaeO  ( j  )  /3 .  -t-etael  (  j  )  ) 

107  pizzi=-l ourby3* ( etaiO ( j ) /3 . +etail ( j ) ) 

108  pizz=pizze+pizzi 

109  qte=-chiprpe(j)*boltz*tel 

110  qti=-chiprpi(j)*boltz'»til 

111  dilux(2, j)=pizz/(den*mi) 

112  dllux(3, j)=qte/(l .Beboltz*den) 

113  dflux(4, j)=qti/(l .5*boltz*den) 

114  dsource(2,j)=-dflux(2,j) ♦ v 1* (denl/den) -drag ( j ) 

115  <l«ource(3, j)=-(boltze(qeilac(j)*(te-ti)+qeOfac(j)*{te-tempO)) 

116  2  +pizze*vl+qte’*(denl/den))/(1.5*boltz*den) 

117  dsource(4, j)=-(boltz*{qeiiac(j)e(ti-te)+ci01ac(j)*(ti-temp0)) 

118  2  +pizzi*vl+qti*(denl/den))/(1.6*boltz*den) 


7.5  Results  of  Afodifl  cat  to  ns  to  Terms  in  the  Equations 

With  the  term  /?’"  added  into  the  expression  for  S^,  the  model  ran  for  all  25  cycles  and 
converged  to  a  solution.  For  the  test  case  parameters  there  was  no  noticeable  difference  between 
the  results  with  the  term  /2‘"  and  without  the  term  R*’*. 
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When  run  with  the  extra  temperature  gradient  terms  removed  from  the  differential  equations 
for  temperature,  the  model  evidently  suffered  a  strong  interaction  with  one  of  the  outflow  bound¬ 
ary  conditions  which  led  to  numerical  instabilities.  The  symptoms  were  an  increase  in  electron 
temperature  of  over  two  orders  of  magnitude  in  the  40**’  (rightmost)  integration  bin  which  led  to 
a  decrease  in  the  integration  time  step  in  dirk2  beyond  its  lower  bound  of  1.0  x  10“**  sec. 

Running  the  model  with  the  extra  velocity  gradient  term  removed  caused  rather  different 
behavior.  The  model  ran  for  over  16|  hours  but  completed  only  5  time  cycles.  At  the  onset  of 
numerical  instability  during  the  fifth  time  cycle,  there  was  a  sharp  discontinuity  in  the  density  profile 
and  several  smaller  discontinuities  in  the  velocity  profile.  Evidently,  these  forced  the  integrator  to 
smaller  and  smaller  time  steps  until  it  exceeded  its  lower  bound  of  1.0  x  10“**  sec. 

7.6  Conclusions  from  the  Results  of  Modifications 

Since  there  is  no  noticeable  difference  between  runs  that  include  and  runs  that  exclude  the 
term  R'",  it  is  apparently  negligible  for  the  parameters  of  the  published  results. 

The  sharp  increase  in  electron  temperature  in  the  outermost  bin  observed  during  the  lun 
without  the  extra  temperature  gradient  terms  was  typical  of  the  onset  of  numerical  instabilities 
in  the  model.  Future  work  with  this  model  should  isolate  the  boundary  condition  responsible  and 
investigate  the  nature  of  the  interaction. 

The  model  contains  user-input  stability  factors  in  the  terms  for  the  dissipative  fluxes.  It  may 
be  possible  to  adjust  the  stability  terms  built  into  the  model  to  make  the  modified  equations  stable. 
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Mil.  Enhancements  to  the  Model. 


To  enhance  the  information  provided  by  the  model  and  to  allow  comparisons  (1)  between  this 
model  and  experiment  and  (2)  between  this  model  and  other  models,  the  output  weis  extended  to 
calculate  four  additional  variables;  the  electron  flux,  the  electric  field,  the  electrical  potential,  and 
the  electron  pressure. 

8.1  Eltciron  Flux 

In  addb  ion  to  the  four  basic  macroscopic  quantities  that  the  model  calculates,  the  pl^tsma 
density,  drift  velocity,  electron  temperature,  and  ion  temperature,  another  basic  property  of  a 
plasma  is  the  electron  flux,  F*  ,  defined  by  F*  =  [21:83].  The  model  assumes  (1)  a  quasineutral 

plasma  so  that  n*  =  n‘  and  (2)  ambipolar  diffusion  so  that  V^-'  =  Vi.  Thus  Ff  =  n*  Vj*  =  n'Vi  —  FJ, 
so  the  electron  and  ion  fluxes  are  equal.  The  continuity  equation.  Equation  46  in  Appendix  A,  shows 
that 

^  _ 

di  dxf  ~ 

Because  for  every  electron  created,  there  is  an  ion  created,  5'^,,  =  In  steady  state,  there  is 

no  change  in  time,  so  in  one  dimension: 


Thus  the  slope  of  the  electron  flux  curve  givf’s  a  measure  of  the  net  production  in  the  plasma. 

Figure  19  shows  the  particle  flux  for  the  steady  state.  The  slope,  and  hence  the  net  production 
Scoii.  is  roughly  constant  from  z  =  0  to  z  =  20.4  cm.  From  j  =  20.4  cm  to  r  =  25cm,  the  slope  is 
shallower,  corresponding  to  a  somewhat  smaller  value  of  Scoii-  Because  in  a  real  tandem  volume 
ion  source  the  magnetic  filter  field  confines  the  ionizing  primary  electrons  to  the  production  region, 
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Particle  Flux 


Position*  z  (cn) 

Figure  19.  Charged  Particle  Flux. 

one  would  intuitively  expect  the  net  production  to  be  smaller  in  the  extraction  region.  In  this  case 
the  results  from  the  model  behave  as  intuiiion  suggests. 

However,  as  Section  6.3  demonstrated,  the  effect  of  the  magnetic  filter  field  is  not  included 
explicitly  in  the  model,  but  enters  only  indirectly  through  the  transport  coefficients.  This  enhance¬ 
ment  to  the  model  does  not  calculate  the  flux  as  the  integral  of  the  production  term  Scoii ,  but  as 
the  product  of  the  density  (in  Figures  11  and  14)  and  the  drift  velocity  (in  Figure  11). 

8.2  The  Electric  Field 

Another  enhancement  to  the  model  calculated  the  electric  field  by  three  different  methods, 
and  graphed  the  results  for  comparison.  The  first  method  solved  for  the  electric  field  from  the  ion 
momentum  equation,  the  second  method  solved  for  the  electric  field  from  the  electron  momentum 
equation,  and  the  third  method  solved  for  the  electric  field  from  the  difference  of  the  two  momentum 
equations. 
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The  electric  field  in  terms  of  the  ion  equation  [35]  is: 


enEj  =  m  nj  V^)  -  H - 


dV,\  .  .  dV,\ 


HJ+i)  - 


where  subscript  j  indicates  a  quantity  calculated  in  the  integration  bin,  and  the  numbers  of 
the  terms  correspond  to  the  numbers  of  the  terms  in  Equation  9  in  Chapter  VI. 

Solving  for  the  electric  field  from  the  electron  momentum  equation  yields: 


enEj  =  -m‘nj  V,)j  -  - 

"  -  S 
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Similarly,  the  electric  field  in  terms  of  the  difference  of  the  two  momentum  equations  is: 


enEj  =  ^in^’ -ni")nj  Vz)j 
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Figure  20  gives  the  results  of  the  three  different  methods  of  calculating  the  electric  field  and  an 
additional  graph  of  the  average  of  the  three  methods.  The  general  shape  of  each  of  the  fou'  graphs 
is  similar:  a  slow  linear  rise  in  the  production  region,  a  sharp  and  abrupt  peak  in  the  magnetic  filter 
region,  where  in  a  real  source  charged  particles  are  trapped  or  deflected  by  the  magnetic  field,  and 
a  return  to  the  slow  linear  rise  in  the  extraction  region.  The  lower  left  graph  in  Figure  20,  labelled 
“(ion),”  shows  the  electric  field  derived  from  the  ion  equation.  It  is  smooth  and  regular  and  is 
dominated  by  the  ion  drag  term  (Term  8)  [35].  In  contrast,  the  lower  right  graph  in  Figure  20, 
labelled  “(electron)”  and  derived  from  the  electron  momentum  equation,  is  jagged  and  irregular. 
It  is  dominated  by  the  pressure  term  (Term  3).  Why  there  should  be  such  a  large  difference  in 
the  appearance  of  the  graphs  is  unknown,  as  is  the  reason  for  the  sharp  downward  spike  at  20.3 
cm.  The  upper  right  graph  in  Figure  20,  labelled  “(difference),”  shows  the  electric  field  calculated 
from  the  difference  of  the  two  momentum  equations.  The  appearance  of  this  graph  is  intermediate 
between  the  electron  and  ion  graphs.  The  upper  left  graph,  labelled  “(average),”  is  an  average  of 
the  three  methods  of  calculating  the  electric  field. 
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E-Field  (difference) 
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Figure  20.  The  Electric  Field.  Counter-clockwise  from  lower  left:  the  elctric  field  calculated 
from:  (1)  the  ion  momentum  equation  (“ion”),  (2)  the  electron  momentum  equation 
(“electron”),  (3)  the  difference  of  the  two  momentum  equations  (“difference”),  and  (4) 
an  average  of  the  three  methods  (“average”).  The  filament  cathodes  are  located  at 
2  =  0,  and  vertical  lines  mark  the  1/c  points  of  the  magnetic  filter  field. 
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8.S  The  Poieniial 


Since  the  relationship  between  the  electric  field  in  the  direction  Ei  and  the  potential  $ 


is 


Ei  = 


dxi 


the  potential  in  the  one-dimensional  model  is 


Therefore,  summing  the  contributions  of  the  electric  field  from  each  bin  from  j  =  0  to  j  =  jo 
approximates  the  potential  at  the  location  of  bin  jo-  Figure  21  compares  the  potential  calculated 
from  the  one-dimensional  model  with  a  graph  of  experimental  values  of  the  plasma  potential  from 
Johnson  [18:168].  The  enhancement  to  the  model  calculates  the  electrical  potential  from  each  of 
the  four  methods  of  calculating  the  electric  field.  Figure  21  displays  the  potential  calculated  from 
the  electric  field  derived  from  the  ion  momentum  equation  (labelled  “ion”  in  Figure  20)  because 
the  graph  of  that  method  of  calculating  the  electric  field  is  the  smoothest.  Because  the  process  of 
integration  itself  is  a  smoothing  process,  the  potential  curves  calculated  from  the  other  methods 
are  indistinguishable  from  the  one  displayed  in  Figure  21. 

The  model’s  calculation  of  the  potential  appears  to  agree  well  both  qualitatively  and  quan¬ 
titatively  with  the  trend  of  experimental  data.  Where  the  model  calculates  a  l^volt  in  potential 
from  the  cathode  to  the  extraction  anode,  experiment  shows  a  drop  of  just  under  2  volts.  Since  the 
potential  is  the  integral  of  the  electric  field,  the  zero  point  depends  on  the  choice  of  the  constant 
of  integration.  Thus  the  potential  is  a  relative  quantity,  and  the  scale  values  are  arbitrary.  For 
simplicity,  this  enhancement  chose  an  integration  constant  of  zero. 
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Figure  21.  Plasma  Potential:  Experimental  Data  (Top)  from  Johnson  [18:168]  Compared  to  Cal¬ 
culation  by  the  Model  (Bottom).  In  both  graphs,  the  filament  cathodes  are  located  at 
r  =  0,  and  vertical  lines  mark  the  1/c  points  of  the  magnetic  filter  field. 
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8.4 


The  Electron  Pressure 


Another  important  parameter  in  an  ion  source  is  the  electron  pressure,  p*  ,  which  is  defined 
by 

p*  =  n*r* 

Johnson  [18:163]  explains  why  this  is  an  important  parameter: 

If  small  numerical  constants  from  the  implied  average  over  the  electron  distribution 
function  are  ignored,  this  parameter  can  be  said  to  represent  the  energy  density,  in 
eV  cm~^,  of  the  bulk  plasma. 

Figure  22  shows  results  of  calculations  of  the  electron  pressure  from  the  one-dimensional  model 
compared  with  experimental  results  of  the  electron  pressure  from  Johnson  [18:164]. 

The  general  trend  from  high  electron  pressure  in  the  production  region  to  low  electron  pressure 
in  the  extraction  region  that  the  model  calculates  agrees  with  the  trend  in  the  experimental  data; 
however,  the  shape  and  magnitude  of  the  decline  are  different® .  The  magnitude  of  the  decline 
that  the  model  calculates  is  small  (just  over  15%),  where  the  magnitude  of  the  decline  in  the 
experimental  data  is  almost  two  orders  of  magnitude. 


’  Coincident2dly,  the  shape  of  the  electron  pressure  curve  that  the  model  calculates  is  very  similar  to  the  logarithm 
of  the  experimental  data. 
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Electron  Pressure  (Experimental) 
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Figure  22.  Electron  Pressure  (Energy  Density):  E.xperimental  Data  (Top)  from  Johnson  [18:164] 
vs.  Calculations  by  the  Model  (Bottom).  In  both  graphs,  the  filament  cathodes  are 
located  at  r  =  0,  and  vertical  lines  mark  the  1/e  points  of  the  magnetic  filter  field. 
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IX.  Conclusions  and  Recommendations 


9.1  Conclusions 

The  one- dimensional  plasma  fluid  model  model  studied  in  this  thesis  hcis  a  lot  of  potential: 
The  ability  of  the  method  of  moving  finite  elements  to  concentrate  computational  power  where  it 
is  needed  is  attractive.  The  derivation  of  most  of  the  equations  upon  which  the  model  is  based 
hcis  been  verified.  The  FORTRAN  code  is  well-written  and  portable,  and  it  is  designed  to  allow 
additional  reactions  to  be  included  and  additional  parameters  to  be  calculated  with  only  minor 
changes.  On  the  other  hand,  the  utility  of  the  model  is  limited  by  several  factors. 

(1)  The  results  out  of  the  current  model  are  unphysical.  The  model  predicts  plasma  den¬ 
sity  higher  in  the  extraction  region  than  in  the  production  region,  opposite  from  experimental 
results.  Moreover,  it  predicts  an  electron  temperature  almost  an  order  of  magnitude  larger  than 
that  experimentally  observed. 

(2)  The  simplification  of  the  reaction  chemistry  used  to  develop  the  model  causes  the  model  to 
neglect  the  principal  paths  of  ionization  and  recombination.  The  resulting  pleisma  balance  between 
ionization  and  outflow  is  incorrect,  where  the  correct  balance  would  be  between  ionization  and 
recombination. 

(3)  The  model  neglects  wall  losses. 

(4)  There  are  inconsistencies  in  the  derivations  of  some  equations,  leading  to  extra  terms  in 
the  momentum  and  temperature  equations. 

(5)  The  method  used  to  calculate  the  energy  flux  at  the  left  boundary  makes  the  final  solution 
depend  on  the  choice  of  the  initial  guess  at  a  solution. 

(6)  For  an  unknown  reason,  multiple  runs  of  the  model  do  not  reproduce  the  same  results. 

(7)  When  terms  in  the  equations  are  modified  or  when  initial  conditions  are  varied  outside  of 
a  relatively  narrow  range,  the  model  develops  numerical  instabilities. 


If  these  limitations  can  be  remedied,  the  code  could  become  a  very  useful  tool  in  the  study 
of  volume  tandem  magnetic  multicusp  ion  sources. 

Although  the  reason  that  the  model  predicts  a  plasma  density  that  trends  opposite  to  experi¬ 
ment  is  still  unknown,  the  incorrect  electron  temperature  can  be  explained  in  terms  of  the  processes 
occuring  within  the  model.  The  model  is  apparently  driven  by  an  outflow  boundary  condition  at 
the  extraction  end.  Because  the  model  limits  the  production  of  plasma  to  the  ionization  of  atomic 
hydrogen  H ,  and  because  it  neglects  the  distinction  between  hot  primary  electrons  and  thermalized 
secondary  electrons,  the  outflow  boundary  condition  requires  a  greater  production  than  secondary 
electrons  can  achieve  at  experimentally  observed  temperatures  of  about  leV.  The  model  compen¬ 
sates  by  raising  its  electron  temperature  to  about  7  eV'  to  achieve  the  required  production.  The 
reaction  rate  for  the  ionization  of  atomic  hydrogen  rises  almost  six  orders  of  magnitude  between 
1  eV  and  7  eV,  and  the  production  matches  the  required  outflow  at  this  point. 

The  population  of  hot  primary  electrons  accomplishes  most  of  the  ionization  in  a  real  ion 
source  through  ionization  reactions  with  both  H  and  H2-  Moreover,  the  loss  mechanism  in  a  real 
source  includes  recombination  reactions  with  and  as  well  bs  H"*  .  The  recombination  rates 
for  the  dissociative  recombination  reactions  with  H2  and  are  six  orders  of  magnitude  larger 
than  the  radiative  recombination  rate  with  so  that  the  balance  in  a  real  source  is  between 
production  and  loss.  By  neglecting  the  principal  plasma  loss  mechanisms,  the  model  achieves  an 
incorrect  balance  between  production  and  outflow. 

The  tendency  of  the  model  to  develop  numerical  instabilities  made  implementing  modifica¬ 
tions  much  more  difficult  than  anticipated.  This  tendency  slowed  the  investigation  and  limited  the 
achievement  of  the  objectives  outlined  at  the  beginning.  However,  further  study  may  provide  the 
needed  breakthrough  to  overcome  the  limitations  that  this  study  has  identified. 
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9.2  Recommendations 


A  number  of  recommendations  for  further  study  follow  directly  from  the  conclusions  above. 
Although  the  model  does  show  a  great  deal  of  promise,  it  has  many  limitations,  and  any  future 
investigator  should  weigh  carefully  the  costs  and  benefits  of  continued  investigation. 

9.2.]  Reproducible  Results.  The  most  important  avenue  for  further  investigation  into  the 
model  is  to  determine  the  cause  of  the  variation  between  different  runs  of  the  model  with  identical 
inputs.  Until  the  cause  of  this  behavior  can  be  identified  and  corrected,  the  model  can  be  of  little 
use. 


9.2.2  Boundary  Conditions.  The  boundary  conditions  seem  to  be  one  of  the  main  keys 
to  the  performance  of  the  model.  Further  research  should  more  explicitly  identify  the  boundary 
conditions  and  where  they  are  applied  in  the  code.  Identification  of  the  boundary  conditions  and 
how  they  affect  the  code  would  allow  a  greater  variety  of  modifications,  which  would  hopefully  yield 
insight  into  the  physics  behind  the  model. 

9.2.3  Stability  Terms.  Another  important  area  for  further  investigation  is  the  effect  of  the 
stability  terms  introduced  to  accelerate  convergence  of  the  fluid  equations.  Understanding  the 
behavior  of  these  terms  may  provide  the  means  to  overcome  the  tendency  of  the  model  to  develop 
numerical  instabilities  when  the  terms  of  the  equations  are  modified  and  when  the  initial  conditions 
are  varied  outside  of  a  narrow  range. 

9.2.4  Wail  Losses.  An  accurate  model  of  a  tandem  volume  ion  source  must  account  for 
wall  losses.  It  may  be  possible  to  include  a  wall  loss  term  developed  in  a  manner  similar  to  that  in 
Johnson  [18:41-42]. 

9.2.5  Energy  Input.  A  relatively  simple  enhancement  to  the  model  would  be  a  modified 
energy  input  term  based  on  the  discharge  parameters  of  current  I  and  voltage  V.  Section  5.2 
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pointed  out  that  the  energy  input  should  be  into  the  entire  production  volume  V,  not  just  into  the 
first  integration  bin,  because  the  length  of  the  electron  mean  free  path  is  of  the  order  of  the  length 
of  the  device.  Thus  the  energy  input  term  would  have  the  form: 


Cin  oc 


IV 

V 


9.2.6  Inclusion  of  a  B-Field  Term.  In  the  momentum  transport  equation,  Equation  1.14  in 
Braginskii  [3:209],  there  is  a  term  that  includes  the  magnitude  B  of  the  magnetic  filter  field: 


m“n 


dVf* 

dt 


dXi 


dXi 


+  en 


I  \ 

\  filter  field  term  ) 


+  Ri 


When  Glasser  and  Smith  adapt  this  three-dimensional  equation  into  their  one-dimensional 
Equation  3.2,  the  term  for  the  magnetic  filter  field  disappears.  In  a  one-dimensional  model,  the 
cross  product  of  the  drift  velocity  in  the  z  direction  and  the  magnetic  field  in  the  x  direction  is 
in  the  y  direction  and  thus  out  of  axial  direction  to  which  the  model  limits  itself.  So  the  only 
remaining  effect  of  the  magnetic  filter  enters  through  the  transport  coefficients.  However,  there  are 
velocities  in  the  perpendicular  directions  that  cannot  be  neglected  (charged  particles  orbit  when 
they  reach  the  magnetic  filter).  It  seems  that  the  code  should  include  the  effect  of  the  filter  field 
explicitly.  It  should  be  possible  to  include  its  effect  by  following  a  diffusional  derivation  like  that 
in  Chapters  6,  7,  and  9  of  Golant  [10]. 
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Appendix  A.  Derivation  of  the  Moment  Equations 


This  appendix  is  based  on  a  course  handout  which  summarized  a  similar  derivation  in  Chap¬ 
ter  6  of  Holt  and  Haskell  [16:152-176]. 


A.l  Introduction 

By  analogy  to  moments  of  force  acting  on  a  lever  arm,  we  obtain  velocity  moments  of  the 
Boltzmann  equation  by  multiplying  the  Boltzmann  equation  by  the  n'*  power  of  the  velocity  and 
then  integrating  over  the  velocity  space.  The  powers  of  the  velocity  generally  chosen  are  v° ,  mv^, 
and  to  correspond  to  the  respective  physically  measurable  quantities  of  number,  momentum, 

and  energy. 


A. 2  Definitions 

In  order  to  derive  the  velocity  moments  c  the  Boltzmann  equation,  we  first  need  some  basic 
definitions.  The  expectation  value  {<j>)  of  some  quantity  <i>  is  defined  as: 


{<!>)  = 


/  <l>n°f^dv°  _  J_ 


(18) 


where  v  is  the  vector  velocity  of  species  a,  /“(v,x,t)  is  the  distribution  function  for  a  single 
particle  of  species  a,  n“(x,  t)  =  /  n®/“(v,x,<)dv“  is  the  density  of  particles  of  species  q  at  point 
X  and  time  t  [21:83],  and  n®  =  [21:79,  316]  is  the  average  density — the  ratio  of  the  total 

number  of  particles  iV®  of  species  a  to  the  total  volume  V  occupied  by  those  particles.  From  this 
definition  of  the  expectation  value  it  follows  that,  if  a  and  6  are  constants  and  g  and  h  are  functions 
of  V,  then 


(a)  =  a 

(19) 

(ag)  =  a{g) 

(20) 
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(ag  +  bh)  =  (ag)  +  (bh)  =  a{g)  +  b{h) 


(21) 


The  telocity  v  can  be  decomposed  into  two  components,  the  mean  Jfelocity  V  and  the  pecj/liar 
velocity  u,  so  that  for  each  species  a, 

<  =  +  ut  (22) 

where  i  is  an  index  defined  so  that  i  =  x,  y,  z.  Since  the  peculiar  velocity  is  random,  it  averages 
to  zero: 

(ti.)  =  0  (23) 

Because  T(x,<)  is  noi  a  function  of  velocity  u,  the  average  velocity  V  can  come  outside  of  ’  e 
integral  in  Equation  18  so  that 

(Vi)  =  VS  (24) 

We  can  apply  the  results  in  Equations  23  and  24  to  obtain  the  expectation  value  of  the  product 

UiVj-. 

(uiVj)  =  Vj{ui)  =  0  (25) 

Moreover,  from  Equations  22,  21,  23,  and  24,  (v,)  =  (Vi  +  u,)  =  (V^)  +  (u;)  =  VS  so  that 

(vi)  =  V-,  (26) 


which  is  just  the  definition  of  the  average  velocity. 

In  order  to  derive  the  moments  of  the  velocity  distribution,  we  also  need  the  relationships 
(viVj)  and  (viVjVj).  If  we  break  the  velocity  into  its  components  in  the  manner  of  Equation  22,  the 
first  product  becomes 

(viVj)  =  ((v;  +  Ui)(Vj  +  Uj))  =  {(ViVj  +  UjVi  +  UiVj  +  u,U;)) 
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Since  the  expectation  value  of  a  sum  is  just  the  sum  of  the  expectation  values, 


{viVj)  =  {ViVj)  +  (ujVi)  +  {uiVj)  +  {uiUj) 


Applying  Equation  24,  we  obtain 


(viVj)  =  ViVj  +  Vi(uj)  +  Vj{ui)+{uiUj) 
=0  =0 


Recognizing  that  the  middle  two  terms  are  zero  from  Equation  25,  our  final  relation  becomes 


(viVj)  =  ViVj  +  {uiUj) 


(27) 


The  triple  product  can  be  analyzed  in  the  same  way.  Applying  Equation  22  to  break  the 
velocity  into  its  components  gives 


{viVjVj)  =  ((V^-  +  u,)(v;'  +  +  Ui))  = 


{ViVjVj  +  ViVjUj  +  ViUjVj  +  ViUjUj  +  UiVjVj  +  UiVjUj  +  UiUjVj  +  UiUjUj) 


Applying  Equations  21  and  24  and  collecting  terms  yields 


{viVjVj)  =  ViVjVj  +  2ViVj  (uj)  +VjVj  («.}  -^2Vj{uiUj)  +  Vi{xtjUj)  +  (u<«,Uj) 

=0  =0 

Equation  23  shows  us  that  once  again  two  of  the  terms  are  zero,  so  our  final  expression  for  the 
expectation  value  of  the  triple  product  is 


(viVjVj)  =  ViVjVj  +  2Vj{uiUj)  +  Vi{ujUj)  +  {uiUjUj) 


(28) 
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A. 3  The  Generalized  Equation  of  Change 

If  /"{v,x,I)  is  the  distribution  function  of  a  single  particle  of  species  a,  and  n®  is  the  average 
density  of  particles  of  species  a,  as  defined  above,  then  n"/“(v,x,<)  is  the  density  distribution  of 
particles  of  species  a  with  velocities  in  the  range  v  to  v  +  rfv  and  positions  in  the  range  x  to  x  +  dx. 
The  Boltzmann  equation  for  particles  of  species  a  has  the  form: 


where 


«.  =  ^  =  ^(Ei  +  kujkvfBt])  =  — (E  +  i[v„ 


Sf^\ 

is  the  acceleration  due  to  external  forces  acting  on  the  particles,  )  is  the  time  rate  of  change 

/  coll 

due  to  collisions,  Z  is  the  charge  state  (defined  as  Z  =  —1  for  electrons  and  Z  =  +1  for  hydrogen 
ions),  and  e.jt  is  the  alternating  unit  tensor,  defined  by  [16:36] 


1  if  the  indices  are  cyclic:  xyz,  yzx,  zxy. 

€ijk  =  '  —1  if  the  indices  are  anticyclic:  zyx,  yxz,  xzy. 

0  if  any  two  or  more  indices  are  equal . 


Let  ^  be  a  generalized  power  of  the  velocity,  so  that  we  can  later  set  <f>  =  1,  <l>  =  mv,  and 
d)  =  ^mv^.  Multiplying  the  Boltzmann  equation.  Equation  29,  by  <t>  and  integrating  over  the 
velocity  space,  we  get 
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which  uses  the  standard  shorthand  notation  for  integration  over  the  velocity  space: 


j  id\  =  J  J  J  ^  ~  J  ^ 


t  =  x,y,z 


(33) 


for  some  arbitrary  quantity  ^  . 

To  develop  the  generalized  equation  of  change  in  terms  of  (^),  we  apply  the  chain  rule  and 
rearrange  terms.  In  the  case  of  the  first  term,  the  chain  rule  gives 


sr--"  f  IT 


(34) 


where  the  term  corresponding  to  the  first  term  in  the  integrated  Boltzmann  equation  is  indicated 
with  an  underbrace. 


Solving  for  the  term  labelled  “1”  on  the  right  hand  side  of  Equation  34,  we  can  rewrite  the 
first  term  of  the  integrated  Boltzmann  equation.  Equation  32,  as 


/ 


dt 


1 


di 

la 


Because  we  are  looking  at  the  temporal  evolution  of  the  system  in  terms  of  a  phase  space  of 
&N  dimensions,  where  N  is  the  number  of  particles  in  the  system,  all  of  the  primary  coordinates 
(x,  y,  z,  Vx,  Vy,  Vx,  and  t)  are  independent  from  each  other.  This  independence  allows  us  to  bring 
the  time  derivative  outside  of  the  integral  in  Term  la  to  obtain 


1 


la 


U 
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From  the  definition  of  the  expectation  value  ^  in  Equation  18  above,  we  get  the  relation; 


dip  . 


By  substituting  into  the  above  relation  for  4>t  we  obtain  another  relation: 
at 


Applying  these  two  relations  to  the  expression  for  the  first  term  of  the  integrated  Boltzmann 
equation  gives 

(35) 

1  la  lb 

where  the  integrals  really  haven’t  been  solved,  but  have  merely  been  hidden  by  the  bracket  notation 
for  expectation  values. 

The  second  term  of  the  integrated  Boltzmann  equation  lends  itself  to  similar  analysis.  In  this 
case  the  chain  rule  gives 


dxf 


(36) 


Solving  for  the  term  marked  “2”  and  using  the  resulting  expression  to  expand  the  second  term  of 
the  integrated  Boltzmann  equation  gives 

Ji.n-vf^dy-  =  I  =  j  f 


2a 


2b 
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Once  again  we  can  apply  the  independence  of  the  v,  and  x,-  coordinates  to  bring  the  partial  differ¬ 
entiation  outside  the  integral  in  Term  2a: 


Substituting  into  the  definition  of  the  expectation  value  in  Equation  18  gives  a  relation  for  {vf<l>) 


and  another  for  (vf 


We  can  rewrite  the  second  term  of  the  Boltzmann  equation  in  terms  of  these  two  relations  to  get 


the  final  form  of  that  term: 


where  Term  2a  is  analogous  to  Term  la,  and  Term  2b  to  Term  lb. 

To  derive  a  term  like  the  third  term  in  the  integrated  Boltzmann  equation,  Equation  32,  we 
can  apply  the  chain  rule  to  the  following  expression: 


—(4>n  aj  )-n  a,/  ^  a, — 


87 


Solving  for  the  term  labelled  “3”  in  Equation  38  and  using  the  resulting  expression  to  expand  the 
third  term  of  the  integrated  Boltzmann  equation  gives 


^ .  V  ^ 

3 


=  /  ^(<in»afr)dv“  + 

> - 1 - ' 

3a 


/ 


d<^ 


36 


(39) 


Looking  at  Term  3a,  we  see  that 


I 


af4>n°f^ 


1+00 
1  —  00 


Because,  as  Vi  goes  to  infinity,  the  distribution  function  goes  to  zero  more  strongly  than  <^  =  1, 
4>  =  mv,  or  <^  =  imv®  go  to  infinity,  the  value  of  term  3a  is  zero. 

Likewise,  we  can  see  that  the  value  of  term  3c  is  zero  by  evaluating  the  partial  derivative 
5a? 

Equation  30  gives  a  value  for  the  acceleration  which  we  can  substitute  into  the  derivative  to 
separate  it  into  electric-field  and  magnetic-field  components: 


_d 

dv 


Ze  dEi 


m®  avf 

E— field  component 


Ze  d 
m°c  dvf 


B— field  component 


First,  consider  the  electric-field  component.  Since  £’,  is  not  a  function  of  v,. 


dEi 

dvf 


=  0 


so  the  value  of  the  electric-field  component  is  zero. 
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Next,  consider  the  magnetic-field  component.  Since  neither  Bt  nor  fjjj,-  are  functions  of  r,, 


a 

dvf 


UVJ 

[cijkvfBt]  = 


Looking  at  the  partial  term  on  the  right,  we  see  that 


dvf  1  i  =  j 

0  i  j 


and  there  is  no  contribution  when  i  j.  But  from  the  definition  of  the  alternating  unit  tensor  f.-jt 

in  Equation  31,  when  i  =  j,  then  =  0.  So  when  i  =  j,  the  magnetic-field  component 

dvf 

IS  zero  because  e,;*  =  0,  and  when  i  ^  j,  the  magnetic  field  component  is  zero  because  — —  =  0 

dvf 

Therefore,  both  the  components  of  Term  3c  are  zero,  so  the  value  of  Term  3c  is  also  zero. 

Only  Term  3b  remains  in  the  third  term  of  the  integrated  Boltzmann  equation.  The  definition 
of  the  expectation  value  allows  us  to  rewrite  this  term  more  compactly: 


36 


Substituting  the  values  derived  for  Terms  3a,  3b,  and  3c  back  into  Equation  39  gives  the  final 
expression  for  the  third  term  of  the  integrated  Boltzmann  equation: 


dvf 


(40) 
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where  Term  3b  is  of  the  same  form  as  Terms  lb  and  2b. 
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Substituting  the  expressions  in  Equations  35,  37,  and  40  into  Equation  32,  the  integrated 


Boltzmann  equation  becomes: 


1 


_d 

dt 


If  we  collect  together  the  “a”  terms  and  the  “b”  terms  in  Equation  41,  we  get  what  is  known  as 
the  generalized  equation  of  change: 


la  2a 

/  16  26  36  \ 

4 

=  } 

- ^ 

^  (42) 

/  coll 

1  2 

V  3  '  'k  '  '  s  ') 

6 

where  the  overbraces  identify  the  terms  of  the  integrated  Boltzmann  equation,  and  the  underbraces 
identify  the  terms  of  the  generalized  equation  of  change. 

Holt  and  Haskell  [16-.154-155]  identify  the  physical  significance  of  the  terms  of  the  generalized 
equation  of  change:  Term  1  is  the  “local  time  variation  of  the  mean  value”  of  4>  [16:154].  Term  2 
represents  the  part  of  that  local  time  variation  caused  by  “molecules  streaming  in  and  out  of  the 
volume  element”  [16:154].  Terms  3,  4,  and  5  represent  the  part  of  the  local  time  variation  due  to 
the  dependence  of  <i>  on  “the  time,  the  position  of  the  particles,  and  the  external  forces”  [16:155], 
respectively.  Finally,  Term  6  represents  the  part  of  the  local  time  variation  caused  by  collisions  with 
particles  of  other  species.  Any  gain  in  the  mean  value  of  must  be  reflected  in  a  corresponding 
and  equal  loss  from  the  other  species. 
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A. 4  The  Zeroth  Velocity  Moment  (The  Continuity  Equation) 


To  develop  an  equation  for  number,  we  let  ^  =  1  in  the  generalized  equation  of  change.  Then 

84)  84)  84> 

—  =  0,  =  0,  =  0,  and  (4))  =  1  so  that  the  generalized  equation  of  change  reduces  to 


It  is  customary  to  define  the  term  so  that 

/  -^1  <^v“  =  Sg^ii  =  (Production)  —  (Loss)  (44) 

J  /  coll 

where  the  term  production  refers  to  production  of  charge  carriers,  or  ionization,  and  the  term  loss 
refers  to  the  loss  of  charge  carriers,  or  recombination.  If  we  substitute  from  Equation  44,  and 
recognize  that  Equation  26  shows  that  (n“)  =  VA ,  then  Equation  43  becomes 

+  =  («) 

This  equation  expresses  the  principle  of  conservation  of  number  of  charged  particles  and  is  therefore 
called  the  continuity  equation.  Equation  45  is  the  same  form  used  in  Equation  1.11  of  Braginskii's 
derivation  [3;208],  except  that  Braginskii  eissumes  that  no  particles  are  created  or  destroyed  in 
collisions  so  that  on  the  right  hand  side,  =  0. 

Applying  the  chain  rule  to  the  second  term  on  the  left  side  of  Equation  45  gives 

8xf  *  8xf 
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Rewriting  the  continuity  equation,  Equation  45,  with  the  second  term  expanded  by  means  of  this 
expression  gives: 

+  "  j-®-" 

VVe  can  now  regroup  the  terms  thus: 


dxf 


to  get  another  form  of  the  continuity  equation: 


di" 


(46) 


where 


dt-  at  ’  axf 


(47) 


is  the  definition  of  the  convective  derivative,  ^ 


A. 5  The  First  Velocity  Moment  (The  Momentum  Equation) 

To  obtain  an  equation  for  the  momentum  balance,  we  set  <i>  =  m“v“  and  substitute  into 
Equation  42,  the  generalized  equation  of  change,  to  obtain: 


an“((r7i“rf))  ^  an°(i',°(m"i;/)) 


at 


dxf 


—  n 


/  ^ 

'  M  9x9'  '  ^  •  9„9  ' 
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Since  (<^}  =  7n“(t;“)  =  m^V°  by  Equation  26,  we  can  transform  the  first  term  of  the  gener¬ 
alized  equation  of  change: 

))  =  ^(n°m‘'V°) 

dr  ^  ^  ■>  \ 

1  1' 

From  Equation  20  and  Equation  27  we  can  see  tliat 


=  (vrim^vf))  =  w^ivrvf)  =  m‘^V^Vf  m“(ufu;) 


so  we  can  expand  the  second  term  of  the  generalized  equation  of  change; 


dxf 


2o' 


24' 


Since  the  velocity  apace,  configuration  space,  and  time  coordinates  are  independent,  the  time 
and  configuration  space  derivativej  of  are  z3ro: 


dt 


di 


dxf  dxf 


=  0 


This  means  that  the  third  and  fourth  terms  of  the  generalized  equation  of  change  are  also  zero: 


dvf 


dvf 


If  we  use  the  Kronecker  delta  6ij  with  the  property  that 


1  i  =  j 
0 
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Collecting  the  primed  terms,  we  obtain  a  form  of  the  equation  for  momentum  balance: 


(49) 


where  the  overbraces  indicate  terms  of  the  generalized  equation  of  change,  Equation  42,  and  the 
underbraces  indicate  terms  of  this  form  of  the  momentum  balance  equation. 

Each  of  the  terms  of  the  momentum  balance  equation,  Equation  49,  has  a  physical  interpreta¬ 
tion  [6:80-82].  Term  1  is  the  time  rate  of  change  of  the  momentum  density.  Term  2  is  the  gradient 
of  the  macroscopic  momentum  flow  (momentum  density  times  velocity).  Term  3  is  the  gradient 
of  the  “hidden”  (microscopic)  momentum  flow,  due  to  the  random  thermal  motion  of  molecules. 
Term  4  is  the  change  in  momentum  caused  by  external  forces.  Term  5  is  the  change  in  momentum 
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of  particles  of  species  o  due  to  collisions  with  particles  of  other  species.  Any  gain  or  loss  through 
this  term  will  appear  in  the  momentum  balance  equation  of  another  species. 

It  is  possible  to  rewrite  the  first  two  terms  of  Equation  49  in  terms  of  the  continuity  equation, 
Equation  46.  Using  the  chain  rule  to  expand  each  term  gives 


/ 


dVf  dn° 

n"  m“  -ei-  +  Vf - 

dt  ^  dt 

- - - '  ' - - - 

\  la  n  / 


\  I 
+ 


3V,“  dn'^V° 

n°m°'  +  ni“  1/“ 

'  flr-^  ^  J  dx9 


V 


2a 


26 


Regrouping  the  terms  on  the  right  side  of  the  equation  gives 


dxf 


/Continuity  EquationX 


dV9 

dxf 

dt 

Ta  ^ ) 

where  the  quantity  multiplying  m‘’Vf  on  the  right  hand  side  is  the  left  hand  side  of  the  continuity 
equation,  Equation  45.  If  we  substitute  from  the  continuity  equation,  the  first  two  terms  of  the 
momentum  balance  equation  transform  to: 


V  ■■  ^  ^  “V 


la  &  2a 


16&  2b 


where  the  convective  derivative  —  is  defined  as  in  Equation  47. 


To  derive  an  expression  for  Term  3  in  Equation  49,  we  need  to  look  more  closely  at  the 
kinetic  energy  and  its  relationship  to  pressure.  The  expectation  value  of  the  kinetic  energy  density 
of  particles  of  species  a  and  velocity  vf  may  be  expanded  using  Equation  20  and  Equation  27  to 
obtain: 
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The  first  term  in  this  expansion  may  be  identified  as  the  macroscopic  kinetic  energy,  and  the  second 
term  as  the  microscopic  or  thermal  kinetic  energy. 

From  the  definition  of  the  expectation  value  in  Equation  18,  we  can  write  the  thermal  kinetic 


energy  density  as: 


2  '  '  ’ '  2  / 


Suppose  that  f^{u^)  is  approximately  a  Maxwellian  distribution  [29:492]  [21:5],  so 

/“K)  =  ^exp|— 

/  ni 

where  A  =  rp'j  is  the  normalization  constant  for  the  distribution.  Substituting  this 
Maxwellian  distribution  into  the  expression  for  the  thermal  kinetic  energy  gives; 

1  1  /«n“>4exp{-::!^^}dv“ 

^  .  I-:  - - 

^  ^  /nMexp| — 


Multiplying  the  numerator  by  1  =  —  and  rearranging  terms  enables  us  to  form  the  argument 
of  the  exponent  inside  the  integral  and  bring  the  constants  outside  the  integral,  where  they 


cancel: 


1  TI*^A  J  2T“  exp  S  2’^“  ^  f 

_a  ^  /..a,.a\ _ oo-'a  •*  ^  **  **  I  J 


TTl^U^U^ 

The  change  of  variables  ^  =  — ^a~  implies 


I2m°‘  m°ufuf  ^  „ 
T“  V  27’“ 


This  expression  for  may  be  solved  for  duf  to  give 


/T“  1 
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Since  is  a  constant,  Equation  33  and  Equation  22  show  that 


dv“  =  dvf  =  +  O  =  dV;“  +duf  =  duf 


so  we  can  write 


= 'fS' 


=  0 


di 


Upon  making  the  substitution  of  variables  into  the  integral,  the  remaining  constant  term 


/  T“ 

V  2m“ 


can  come  outside  the  integral,  so  the  expression  for  the  thermal  kinetic  energy  density  becomes: 


=  3n"T" 

2  exp  {-Odv- 

=1 


where  the  factor  of  3  comes  from  summing  over  the  index  %  =  x,y,z.  The  gamma  function  r(C)  is 
defined  by 


and  has  the  property  that 


r(c)=  rt^<-%-*dt  c>o 

Jo 


r(c  + 1)  =  cr(c) 


Recognizing  the  gamma  functions  in  the  numerator  and  the  denominator,  we  can  transform  the 
expression  for  the  microscopic  kinetic  energy  density  to: 


1  irfil  3  3 

n‘‘-m“{ufun  =  =  3n“r“^;^  =  -n“r“  = 


r(^) 


r(t)  2 


which  demonstrates  the  basic  relationship 


(50) 
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between  the  thermal  kinetic  energy  density  and  the  scalar  pressure.  Define  the  quantity  ipfj  as  in 
Equation  1.19  in  Braginskii  [3:210]  so  that 

V'g  =  n“)  (51) 

Then  the  quantity  ipfj  is  a  second-order  tensor  known  as  the  pressure  tensor^ .  Suppose  that  is 
spherically  symmetric,  i.e.,  i  =  j.  Then 

=  n°m°{utul)  =  2n°T°  = 

as  demonstrated  above.  The  spherically  symmetric  part  of  the  pressure  tensor  comprises  its  diagonal 
elements  so  that 

P?  0  0 

^kk  =  0  p“  0 

0  0  p“ 

where  pf  is  the  pressure  in  the  direction.  In  an  isotropic  medium,  p?  =  Py  =  P?  =  p“  where  p“ 
is  the  scalar  pressure.  The  trace  of  V’tt  is  thus 

P?+P?+P“  =  3p“ 

so  that  the  symmetric  part  of  the  pressure  tensor  can  be  expressed  as 

P°^ij  = 

If  the  pressure  tensor  is  not  spherically  symmetric,  as  in  the  ceise  where  a  magnetic  field  is 
present,  it  contains  off-diagonal  terms  and  thus  may  be  written  as  the  sum  of  its  diagonal  and 

'The  symbol  V'  was  not  chosen  to  represent  the  pressure  tensor  because  of  the  abbreviation  p.s.i.,  although  the 
coincidence  is  interesting. 
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ofF-diagonal  elements: 


t/-5=p"^.>  +  n?. 

where  11?  is  the  viscosity  stress  tensor,  defined  by: 

Thus  Term  3  of  Equation  49  can  be  separated  into  two  components,  one  corresponding  to  the  scalar 
pressure,  the  other  to  the  stress  tensor: 

' — 1 - ^ ' 

3 

Equation  30  gives  an  expression  for  the  acceleration  that  we  can  substitute  into  Term  4  of 
Equation  49: 


dxf  dxf  dxf 


4 

Because  Cijk  and  Bk  are  independent  of  Uj,  we  can  bring  them  outside  of  the  expectation  value,  so 
applying  Equation  22  to  break  the  velocity  into  its  components  gives: 


=  Zen^  ^(£7.)  +  6^ 

4 

Since  {tij)  =  0  from  Equation  23,  (V})  =  Vj  from  Equation  24,  and  (£■,)  =  Ei  because  Ei  is 
independent  of  r,-,  the  final  version  of  Term  4  is 
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We  can  use  Equation  22  to  expand  Term  5,  ^4",  into  its  components: 


where  Equation  44  defines  the  net  production,  and 


HI) 

)  eon 


dv" 


(52) 


defines  Rj  ,  the  net  microscopic  momentum  exchanged  in  the  direction  through  collisions  with 
species  a. 

Collecting  the  derived  forms  of  Terms  1  through  5  gives  the  equation: 


If  we  cancel  the  terms  that  appear  on  both  sides  of  the  equation,  we  obtain  the  mo¬ 

mentum  balance  equation.  Equation  1.14  in  Braginskii  [3:209]: 


where  the  overbraces  indicate  the  terms  of  the  original  form  of  the  momentum  balance  equation, 
Equation  49,  and  the  underbraces  indicate  the  terms  of  the  momentum  balance  equation  as  imple¬ 
mented  in  the  model,  Equation  9. 
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A. 6  The  Second  Velocity  Moment  (The  Energy  Equation) 


To  get  an  equation  for  the  energy  balance,  we  begin  with  the  generalized  equation  of  change, 
Equation  42: 


/ 


\ 


d{n-{4>))  a(n°(uf0)) 

at  dxf 

' — ^ 


a  ^^  \  ,  /  a  V 


and  set  4)  equal  to  the  kinetic  energy,  that  is,  4>  =  Because  of  the  independence  of 

coordinates,  the  partial  derivatives  of  <j>  with  respect  to  t  and  x"  are  both  zero: 


d<l> 

It 


1  ^dvfvf 


2"* 


at 


a<f> 

— —  =  -m 

ax?  2 


1  ^avfvf 


ax? 


Therefore,  Term  3  and  Term  4  of  the  generalized  equation  of  change  are  also  both  zero: 


n“O  =  n4m<«(£i|^)  =  0 


'  ax?'  2  '  ’  dx?  ' 


=0 


=0 


The  partial  derivative  of  <t>  with  respect  to  vf  is 


d4>  1 

- =  —  m  -  ■  •' — “ —  = 

dvf  2  avf 


=  < 


» ^  J 


.=i 


avf 


From  the  chain  rule, 


av?  av? 


at;? 


which  is  just  the  familiar  d(t'^)  =  2vdv,  so 


ad>  1  ,dvfvf 
avf  2  avf 


{2vf)  = 
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The  expectation  values  of  <t>  and  vf<i>  ate 


(0)  = 


Substituting  the  values  calculated  above  into  the  generalized  equation  of  change  transforms  it  to 


dt 


1 


2 


5 


6 


where  the  new  quantity  ,  the  exchange  of  energy  due  to  collisions  involving  particles  of  species 
a,  is  defined  by 


Mcon  = 


=  -m  y  Vj  Vj 


6_r\ 

/  coll 


dv“ 


(54) 


We  can  decompose  the  velocity  into  its  macroscopic  and  microscopic  components,  as  in  Equa¬ 
tion  22,  and  substitute  vf  =  -|-  uf ,  for  vf  in  the  equation  of  change.  Using  the  relation  for  the 

product  (vf  )  in  Equation  27,  the  first  term  of  the  transformed  equation  of  change  becomes 


-i- 

dt  \2 


From  Equation  50,  {«“«“)  =  |n“T“,  so  the  first  term  can  be  written  as 


at 


where  =  p°  for  an  ideal  gas. 
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The  second  term  can  be  decomposed  in  a  similar  fashion  using  the  expansion  of  the  product 


{vfv^vf)  in  Equation  28. 


2 


d 

dxf 


2a 


2» 


2c 


2d 


Term  2a  is  the  gradient  of  the  macroscopic  transport  of  kinetic  energy  density  at  the  mean 
velocity  Kj": 

a 

dx? 


/ 

\ 

1 

im®(V®)2 

\ 

m®v;®v/v^® 

_  a 
~  ^ 

n® 

v;® 

2a  / 

1 

Kinetic  Energy 

/ 

Term  2b  represents  the  work  done  by  pressure,  and  Term  2c  represents  the  macroscopic 
transport  of  microscopic  (random  thermal)  energy  at  the  mean  velocity  They  can  be  written 
in  terms  of  the  pressure  tensor  ^<5  =  Substituting  V’i“  into  Term  2b  and  swapping 

indices  in  Term  2c  gives: 


a 

axf 

( 

a 

dxf 

(  \ 

n®m®{u®u®)V,®  +  in®m®  {«;«;>  V;.® 

' - V - '  <2 _ - 

V-“  V,“  +  n®  im®(u“«®)  V;® 

\ 

2c  } 

1 

\  Thermal  Energy  / 

Expanding  t/’i“  =  +  0,®  as  in  the  derivation  of  the  moment  equation  and  using  the  relationship 

5n®m®(u®u®)  =  from  Equation  50,  these  two  terms  become 


=  9^  (^'’•‘‘''''"  +  "5''-“) 
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Contracting  Sij  Vj  over  the  j  index  and  substituting  for  p“  from  the  ideal  gas  law  =  n“T'“  gives 


Term  2d  describes  the  transport  of  microscopic  (random  thermal)  energy  at  the  thermal 
velocity  If  we  define  the  heat  flux  vector  qt  to  be 


qf  =  in“m“(ur«;«p  =  n“(  uf) 

N .  II  I  ^ 

Thermal  Energy 


in  Equation  1.21  in  Braginskii  [3:210],  then  Term  2d  becomes 


Term  5  represents  the  volumetric  power  density  from  the  action  of  external  forces.  Substi¬ 
tuting  for  the  acceleration  from  Equation  30  we  get 


n“(m»a“<)  =  =  Zen-{{Ei  -h  = 

^  ^  m  c 


Since  the  alternating  unit  tensor  f,y*  has  the  property  that  c.yt  =  — for  every  term  [f  i23i'2  , 

where  1,  2,  and  3  represent  some  arbitrary  assignment  of  the  coordinates  x,  y,  and  z,  there  is  a  term 
[f2i3»’f  Bajt)®  =  — 1^123*'“ Sajv®  to  cancel  it.  Therefore,  [(ijk^'j  Bk]vf  =  0,  and  Term  5  becomes 


=  Zen°{EiV°) 
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Since  Ei  is  not  a  function  of  velocity,  it  can  come  outside  of  the  expectation  value,  so  applying 
Equation  26,  we  are  left  with: 

=  Zen°Ei{vf)  =  Zen^- EiV^^  (56) 

' - ^ ' 

5 

Since  Zen"Vj“  is  just  J",  the  macroscopic  electric  current  density  in  the  direction  caused  by 
particles  of  species  o,  Equation  56  is  an  expression  for  the  J  •  E  Joule  heating  of  the  plasma. 

Term  6  represents  the  change  in  energy  density  due  to  collisions  with  particles  of  other 
species.  The  definition  of  in  Equation  54  may  be  expanded  using  the  velocity  decomposition 
v"  =  V'j.®  +  u"  from  Equation  22  to  get 


6<a  66  6c 


Term  6a  contains  as  defined  in  Equation  44,  and  Term  6b  contains  R°  as  defined  in  Equa¬ 
tion  52.  Rewriting  the  expression  for  in  terms  of  these  quantities  gives 

=  IrnOVfVf-Sf^,,  -1-  -h  „ 

where 

defines  the  new  quantity  the  thermal  energy  exchanged  by  the  production  or  loss  of  particles 
of  species  q. 


Collecting  all  of  the  transformations  of  the  terms  of  the  equation  of  change  gives  the  energy 


transport  equation: 


dt 


/, _ 

la 

1  a 

-n 

2 

J  1 

16 


\ 

■  / 

d 

1  j 

dxf 

1 

.  \ 

2a 


2Hi2c 


c+ 


/ 


1 


=  Zen'^EiV^  +  „  +  VfRf  +  „  (58) 


10 


where  the  overbraces  mark  the  terms  of  the  generalized  equation  of  change  and  the  underbraces 
the  terms  of  the  energy  transport  equation.  This  form  of  the  energy  transport  equation  appears  in 
Braginskii  [3:210]  as  Equation  1.20. 
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Appendix  B.  Collision  Processes  in  Hydrogen  Ion  Sources 


This  appendix  contains  a  listing  of  the  primary  reactions  that  can  occur  in  a  hydrogen  plasma, 
grouped  by  the  products  produced  and  the  reactants  destroyed.  In  order  to  standardize  reaction 
numbers  with  those  used  in  the  model,  all  entries,  including  reaction  numbers  K  and  threshold 
energies  Sth.  are  based  on  the  list  in  Smith  and  Glasser  [32]  except  where  cited.  Reactions  are 
listed  in  the  order:  c,  H~ ,  H,  H2,  //'*',  Hf,  and  .  Reactions  with  the  same  first  reactants  are 
listed  in  the  same  order  by  the  second  reaction.  For  more  specific  information  on  reactions  with 
fitted  curves  for  cross  sections  and  reaction  rates,  see  Janev  et  al.  [171. 

B.l  Reaciions  Thai  Release  Electrons 


N 

reaction 

name 

£th 

7 

e  +  H- 

e  ■\  H  +  e 

collisional  detachment 

0.75 

19 

e  +  H 

- 

e  +  H+  +e 

ionization 

13.6 

3 

e  +  ^2  (t)  =  0) 

— 

e+H'  +  H* 

polar  dissociation 

17.32 

24 

e  +  H2 

- 

e  +  H^+e 

ionization 

15.4  [17:52] 

20 

e+H^ 

- 

2e  +  2H+ 

dissociative  ionization 

16.243 

8 

H-  +  H 

- 

e  +  2H 

collisional  detachment  [11:32] 

0.0  [17:230] 

9 

H-  +  H 

^2  (v  =  0)  +  c 

Ekssociative  detachment 

0.0  [17:230] 

58 

H-  +H 

-* 

Hiiv  =  9)  +  c 

(inelastic)  assoc,  detach. 

10 

- 

e  +  H  +  H2 

collisional  detachment 

0.75 

29 

H-  + 

- 

e  +  H^ 

tissociative  detachment 

0,75  [17:222] 

43 

H-  + 

— 

e  +  H* 

associative  detachment 

31 

H  +  H 

— 

e  +  H^ 

associative  ionization  [11:12] 

32 

H  +  H2 

— 

e  +  H* 

associative  ionization  [11:12] 
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B.2  Reactions  That  Capture  Electrons 


K 

reaction 

nume 

— 

£th 

6 

€  +  // 

- 

H~  +  hi/ 

radiative  attachment  [11:43] 

1 

e  +  H2  (t)  =  0) 

- 

H-  +  H 

dissociative  attachment 

3.4 

2 

e  +  H  2(0  =  0) 

- 

H-  +  H' 

dissoc.  attach,  (inelastic) 

13.6 

56 

€  +  II2  (v  =  9) 

H-  +  H 

dissoc.  attach,  (superelasiic) 

0.1  [17:62] 

26 

e  +  H+ 

- 

H  +  hu 

radiative  recombination  [11:47] 

13.6  [17:32] 

4 

e  +  H^(v  =  0) 

H-  + 11+ 

dissociative  attachment 

16 

e  +  lit 

- 

2H 

dissociative  recombination 

0.33 

e+Ht 

H2  +  H 

dissociative  recombination 

0.38 

e  +  Ht 

- 

Ht  +  H- 

dissociative  recombination 
— - - - 1 

B.3  Reactions  That  Produce  H~  Ions 


reaction 

name 

£th 

6 

e  +  H 

-* 

H"  +  hv 

radiative  capture 

1 

e  +  //j  (e  =  0) 

- 

H-  +  H 

dissociative  attachment 

3.4 

2 

e  +  i/2  (v  =  0) 

- 

H-  +H* 

dissoc.  attachment  (inelastic) 

13.6 

3 

e  +  //2  (e  =  0) 

-* 

e  +  H-  +  H+ 

polar  dissociation 

17.32 

56 

c  +  H2{v  =  9) 

- 

H-  +  H 

dissociative  attachment 

4 

e  +  Ht  («^  =  0) 

- 

H-  +  H+ 

dissociative  attachment 

42 

e  +  Ht 

- 

Ht  +  H- 

dissociative  attachment 

22 

H  +  H 

- 

H-  +H+ 

ion  pair  formation  [11:10] 

30 

H  +  H2 

- 

H-  +  Ht 

ion  pair  formation  [11:10] 

21 

H2  +  Ht 
_ 

- 

H-  +  H+  +  Ht 

collisional  dissociation 
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B.4  Reactions  That  Destroy  H~  Ions 


K 

reaction 

name 

£th 

7 

£  +  //- 

- 

e  +  H  +  e 

collisional  detachment 

0.75 

8 

H-  +  H 

- 

e  +  2H 

collisional  detachment  [11:32] 

0.0  [17:230] 

9 

H-  +H 

- 

H‘2  (r  =  0)  +  c 

<issociative  detachment 

0.0  [17:230] 

58 

H-  +H 

- 

Hilv  =  9)  +  e 

associative  detachment 

10 

H-  +  H2 

- 

e^H'hH2 

collisional  detachment 

0.75 

11 

H-  + 

- 

2H 

mutual  neutralization 

29 

H-  A-H+ 

- 

e  +  Ht 

dissociative  detachment 

0.75  [17:222] 

43 

H-  + 

-* 

Ht+e 

associative  detachment 

12 

H-  +  H+ 

- 

H  +  H2 

mutual  neutralization 

B.5  Reactions  That  Broduce  H  Atoms 


H 

reaction 

name 

£th 

■ 

e  +  H- 

e  +  H  A-  e 

collisional  detachment 

0.75 

B 

e  +  Hi 

2H  A- e 

electronic  dissociation 

9.2 

1 

e  +  /fj  (v  =  0) 

- 

H~  A- H 

dissociative  attachment 

3.4 

H 

e  +  i/z  (t'  =  0) 

, 

H-  A- H' 

dissoc.  attachment  (inel^lstic) 

13.6 

56 

e  +  H2{v  =  9) 

- 

H-  A-H 

dissoc.  attachment  (superelastic) 

26 

e  +  H-*- 

- 

HA-hv 

radiative  capture 

13.6  [17:32] 

1 

14 

e  +  H} 

- 

HA-H^+e 

electronic  dissociation 

3.45 

16 

e  +  H+ 

- 

2H 

dissociative  recombination  [11:49] 

0.33 

15 

e  +  H+ 

- 

€A-H  A- Ht 

electronic  dissociation 

6.6 

17 

e  +  H^ 

- 

H2A-H 

dissociative  recombination  [11:49] 

>.38 

8 

H~  A-H 

-- 

eA-2H 

collisional  detachment  [11:32] 

0.0  [17:230] 

10 

H-  +H2 

eA-HA-H2 

(  collisional  detachment 

i  0.75 
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K 


reaction 


name 


£th 


11 

H-+H+  2H 

mutual  neutralization 

12 

H-  +  H^  —  11  +  H2 

mutual  neutralization 

27 

H2  +H+  —  +  H 

charge  transfer  [11:15] 

1.83  [17:142] 

18 

H2  + Hi  -*  H  + 

ion  rearrangement  [11:31] 

0.0  [17:176] 

B.6  Reactions  That  Destroy  H  Atoms 


K 

reaction 

name 

!  ^th 

19 

e  +  H  —  c  +  //+  +  e 

ionization 

13.6 

6 

e  +  H  —  H-+hv 

radiative  capture 

9 

H-  +  H  —  H2{v  =  Q)  +  e 

associative  detachment 

0.0  [17:230] 

58 

H-  +  H  -*  H2{v  =  ^)  +  e 

assoc,  detach,  (inelastic) 

31 

H  +  H  ^  e  +  Ht 

associative  ionization  [11:12] 

22 

H  +  H  —  H-+H+ 

ion  pair  formation  [11:10] 

32 

ff  +  Hi  e  +  H^ 

associative  detachment  [11:33] 

30 

I/  +  H2  -*  D-  +  ff^ 

ion  pair  formation  [11:10] 

33 

H  +  Ht  —  H2  +  Ht 

ion  rearrangement  [11:31] 

B.7  Reactions  That  Eroduce  i/o  Molecules 


K 

reaction 

name 

eth 

40 

e  +  Ht  e  +  H2  +  H-^ 

electronic  dissociation 

17 

e+Ht  -*  H2  +  H 

i 

dissociative  recombination  [11:49] 

0.38 

9 

H-+H  —  ff2(u  =  0)  +  e 

associative  detachment 

0.0  [17:230] 

58 

H-  +  H  —  H2(,v  =  2)  +  e 

associative  detachment  (inelastic) 

12 

H-  +  H+  —  H  +  H2 

mutual  neutralization 

33 

H  +  H}  -  H2  +  H} 

t 

ion  rearrangement  [11:31] 

no 


B.8  Reactions  That  Destroy  H 2  Molecules 


K 

reaction 

name 

£th 

1 

e  +  7/2  (r  =  0) 

- 

H-  +  H 

dissociative  attachment 

3.4 

2 

e  +  H2(v  =  0) 

- 

H-  +  //* 

dissociative  attachment 

13.6 

3 

e  +  H2iv  =  Q) 

- 

e  +  H-  +H-^ 

polar  dissociation 

17.32 

56 

e  +  1/2  (t)  =  9) 

- 

H-  +  H 

dissociative  attachment 

13 

e  +  ^2 

2H  +  e 

electronic  dissociation 

9.2 

24 

e  +  H2 

- 

e  +  +  e 

ionization 

15.4  [17:52] 

30 

H  +  H2 

- 

H-  +  Ht 

ion  pair  formation  [11:10] 

32 

II +  H2 

- 

e-^Ht 

associative  ionization  [11:12] 

27 

H2  +  H+ 

- 

Ht  +  H 

charge  transfer  [11:15] 

1.83  [17:142] 

18 

H2  +  Ht 

- 

H  +  Ht 

ion  rearrangement  [11:31] 

0.0  [17:176] 

21 

H2  +  Ht 

-* 

H-  +H+  + 

collisional  dissociation 

B.9  Reactions  That  Release  Ions 


reaction 

name 

£th 

e  +  H 

- 

e  +  //+  +  c 

ionization 

13.6 

e  +  H} 

- 

2e  +  2H+ 

dissociative  ionization 

16.243 

e  +  Ht 

-* 

H  +  H++e 

electronic  dissociation 

3.45 

4 

e  +  H+  (v  =  0) 

- 

H-  +H+ 

dissociative  attachment 

40 

e  +  H^ 

- 

e+H2  +  H+ 

electronic  dissociation 

22 

H  +  H 

H-  +  H+ 

ion  pair  formation  [11:10] 

21 

H2  + 

- 

H-  +  H+  + 

collisional  dissociation 
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B.IO  Reactions  Thai  Capture  Ions 


K 

reaction 

name 

Eth 

26 

e  +  H+  -* 

H  +  hv 

radiative  capture 

13.6  [17:32] 

29 

H-+H+  — 

e  +  H+ 

associative  detachment 

0.75  [17:222] 

11 

H-+H+  — 

2H 

mutual  neutralization 

27 

H2  +  H+ 

H+  +  H 

charge  transfer  [11:15] 

1.83  [17:142] 

B.ll  Reactions  That  Broduce  Ions 


K 

reaction 

name 

Cth 

24 

e  +  H2 

-* 

e  +  H2  +e 

ionization 

15.4  [17:52] 

15 

e  +  H+ 

- 

e  +  H  +  Ht 

electronic  dissociation 

6.6 

42 

e  +  H^ 

- 

H*  +  H- 

dissociative  attachment 

29 

H-  +  H+ 

- 

e  +  Ht 

associative  detachment 

0.75  [17:222] 

31 

H  +  H 

- 

e  +  Ht 

associative  ionization  [11:12] 

30 

H  +  H2 

- 

H-  +  Ht 

ion  pair  formation  [11:10] 

33 

H  +  H* 

- 

H2  +  H* 

ion  rearrangement  [11  31] 

27 

H2  +  H+ 

- 

Ht  +  H 

charge  transfer  [11:15] 

1.83  [17:142] 

B.12  Reactions  Thai  Destroy  Ions 


K 

reaction 

name 

eth 

20 

e  +  H* 

- 

2e  +  2H* 

dissociative  ionization 

16.243 

14 

e+H* 

- 

H  +  H++e 

electronic  dissociation 

3.45 

4 

e  +  H^(v  =  0) 

- 

H-  +H* 

dissociative  attachment 

16 

e  +  H* 

- 

2H 

dissociative  recombination  [11:49] 

0.33 

12 

H-  +  H* 

— 

H  +  H2 

mutual  neutralization 
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reaction 


name 


B.J3  Reactions  That  Produce  Ht  Ions 


N 

- 1 

reaction 

name 

£th 

43 

H-  +  Ht  -*  H+-i-e 

32 

H  +  H2  —  e  +  H^ 

associative  ionization  [11:12] 

18 

H2  +  H}  —  H  +  Ht 

ion  rearrangement  [11:25] 

0.0  [17:176] 

B.14  Reactions  That  Destroy  H}  Ions 


reaction 

name 

Cth 

15 

c  +  e  +  H  +  Ht 

electronic  dissociation 

6.6 

40 

e  +  H^  —  c  +  Hi  +  H+ 

electronic  dissociation 

42 

e  +  H^  -  Hj+  +  H- 

dissociative  recombination 

17 

e+H^  H2  +  H 

dissociative  recombination  [11:49] 

0.38 

33 

H  +  H}  —  Ho  + 

ion  rearrangement  [11:25] 

B.15  Elastic  Scattering  Reactions 


reaction 


e  +  H  e  +  H  elastic  scattering 

34  e-\-  H2  —*  e  +  II2  elastic  scattering 
50  H~  +  H  — ►  H~  H  elastic  scattering 

5  II~  +  H  — *  H  +  H~  charge  exchange 

35  H~  +  H2  —*  H~  +  H2  elastic  scattering 

25  H~  +  — ►  H'^  H~  double  charge  exchange 


Cth 


0.0 

0.0 

0.0 

0.0  [17:228] 
0.0 
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N 


reaction 


name 


53 

H  +  H 

— 

H  +  H 

47 

H  +  H, 

- 

H  +  H2 

49 

H-¥H+ 

- 

H  +  H+ 

23 

H  +  H-^ 

- 

H-^  +  H 

51 

H  +  H+ 

- 

H  +  H+ 

52 

H  +  H+ 

- 

H  +  Ht 

46 

H2  +  H2 

- 

H2  +  H2 

36 

H2  +  H+ 

- 

H2  +  //+ 

44 

H2  + 

- 

H2  +  Ht 

28 

H2  + 

- 

Ht  +  H2 

45 

H2  +  Ht 

-► 

H2  +  Ht 

eth 

elcistic  scattering  0.0 

elastic  scattering  0.0 

elastic  scattering  0.0 

charge  exchange  0.0  [17:128] 
elastic  scattering  0.0 

elastic  scattering  0.0 

elastic  scattering  0.0 

elastic  scattering  0.0 

elastic  scattering  0.0 

charge  transfer 
elastic  scattering  0.0 


B.16  Ezciiaiion  and  De-exciiaiion  Reactions 


# 

reaction 

name 

Cfh 

37 

e  +  /f3(v  =  0)  -♦  c4-ff2(t>  =  l) 

vibrational  excitation 

0.5  [17:34] 

38 

e  +  (i'  =  0)  — ►  €  +  i/j  (v  =  2) 

vibrational  excitation 

1  [17:36] 

39 

c  4-  //2  (v  =  0)  — *  e  +  ^2  (»'  =  3) 

vibrational  excitation 

1.45 

54 

c  +  /f2  (v  =  0)  — ►  e  +  ^2  (v  =  9) 

vibrational  excitation 

55 

e  +  ^2(r  =  0)  -  H^(,C,B) 

— ►  e  +  H2{v  =  9)  +  hu 

vibrational  excitation 

57 

e  +  H2(v  =  9)  — ►  e  +  H2[v  =  0) 

vibrational  de-excitation 

41 

c  +  ^(ls)  —  e+H(2e) 

electronic  excitation 

10.2  [17:18] 
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Appendix  C.  Glossary  of  Variables 


variable 


Significance 


Arbitrary  or  normalization  constant. 

Momentum  exchanged  due  to  collisions:  A“  = 

B  Magnitude  of  the  magnetic  filter  field. 

Speed  of  light. 

Drag:  d“  =  d  =  — . 

mn 

Electron  charge. 

(Superscript)  Electron:  E.g.,  T‘ ,  F*. 

Electric  field. 

Probability  distribution  function  for  a  single  particle. 

Magnetic  intensity:  yHi  =  B,. 

(Subscript)  First  index:  E.g.,  viVj. 

(Superscript)  Ion:  E.g.,  T*,  F*. 

! 

(Subscript)  Second  index:  E.g.,  ViVj. 

(Subscript)  Integration  bin  number:  E.g.,  7}+i. 

Current  density:  The  current  density  in  the  direction  is  Ji  =  Zen°V°. 
Reaction  rate:  =  ((Thv). 

(Subscript)  Third  index. 

Dimension  of  length. 

Mass:  m  =  m*  +  m'. 

Energy  exchange  due  to  collisions:  dv°. 


variable 

Significance 

n 

Density  n“{x,i)  at  point  x  and  time  t:  n°(x,t)  =  f  n“/"(v,  x,  <)dv“. 

Average  density;  n“  =  N^/V. 

N 

Number  of  particles. 

P 

Pressure:  p“  =  for  an  ideal  gas. 

V 

General  polynomial:  'P{x°). 

Q 

Heat  flux  tensor:  qf  =  ^ti°' m‘‘{ufufuf). 

Q 

Thermal  energy  exchange  due  to  collisions;  f  u°u^h°  11^^°'' 

R 

Momentum  exchanged  due  to  collisions;  /?"  =  f  m“u“n“  ^  dv“. 

s 

Net  production  of  particles  through  collisions;  5°^,,  =  /  n“  dv“ 

t 

Time. 

T 

Temperature  (in  units  of  energy). 

it 

Random  (thermal,  microscopic)  velocity. 

V 

Total  velocity:  Uj  =  V<  +  u,. 

V 

Average  (drift)  velocity. 

V 

Total  volume. 

X 

Argument  to  general  polynomial  P(x°):  x°  =  u°t°'. 

Coordinate  direction  of  magnetic  filter  field. 

V 

Coordinate  direction  perpendicular  to  x  and  z. 

z 

Axial  coordinate  direction. 

Charge  state;  Z  =  —1  for  electrons,  and  Z  =  +1  for  positive  ions. 
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C.2  Greek  and  Hebrew  Letters 


Significance 


First  general  species. 

Second  general  species. 

Particle  flux:  Ff  =  n"V^". 

Gamma  function:  r(C)  =  t^^~^^e~*  dt. 


Kronecker  delta  function  Sij. 


Alternating  unit  tensor  e,jt. 

Energy. 

General  function. 

Coefficient  of  viscosity:  (V{x°)). 

n“T“r'* 

Coefficient  of  thermal  conductivity:  k“  = - (■?(!")). 


1 


Mean  free  path:  Amfn  = 

n®<rK 

Magnetic  susceptibility;  =  B,. 

Collision  frequency:  i/“  = 

Ratio  of  circumference  of  a  circle  to  its  diameter. 

Viscosity  stress  tensor:  =  xpfj  —p°6ij. 

Collision  cross  section. 

Time  between  collisions;  r®  =  — . 

I/O 

Generalized  power  of  the  velocity:  ^  =  1,  ^  =  mv,  and  <!>  =  \mv- 

Electrical  potential. 

Alternative  symbol  for  k. 

Pressure  tensor:  V’ly  =  n®m® {«“«“). 

eB 

Gyrofrequency  of  species  a:  w?  =  - . 

_ m.°C _ 


Unspecified  reaction  number:  E.g., 
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